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WIGNER  DISTRIBUTION  FUNCTION:  RELATION  TO  SHORT-TERM 
SPECTRAL  ESTIMATION,  SMOOTHING,  AND  PERFORMANCE  IN  NOISE 


INTRODUCTION 


The  potential  of  the  Wigner  Distribution  Function  (WDF)  for 
characteri zing  the  short-term  local  time  and  frequency  content  of  a 
transient  waveform  has  been  amply  demonstrated  In  a  series  of  papers;  for 
example,  see  the  recent  publications  [1,2,3]  and  the  extensive  references 
listed  therein.  In  particular,  [1]  contains  numerical  examples  of  the  WDF 
for  rectangularly  gated  linear  frequency  modulation  and  a  version  which  has 
been  smoothed  with  a  square  window  in  the  time-frequency  plane,  in  order  to 
yield  positive  distribution  values.  Here,  we  will  be  concerned  with 
smoothing  so  as  to  minimally  spread  the  WDF,  but  will  not  presume  all  the 
information  that  is  required  for  Implementation  via  [2],  nor  do  we  limit 
consideration  to  a  constant-magnitude  function.  We  will  then  use  the  close 
connection  between  short-term  spectral  estimation  and  smoothed  WDFs  to 
suggest  a  possible  analysis  procedure  and  philosophy  to  extract  Information 
about  a  given  waveform  without  an  extensive  search  in  waveform  parameters. 
Finally,  the  performance  of  a  particular  WDF  estimator  in  the  presence  of 
additive  noise  will  be  analyzed,  both  in  terms  of  bias  and  variance. 

This  report  summarizes  and  compiles  many  of  the  results  in  the 
publications  noted  above  in  a  unified  framework  and  notation.  Also, 
numerous  examples  are  presented  in  the  various  sections  of  this  report  to 
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illustrate  and  bring  out  some  of  the  fundamental  concepts  and  limitations  of 
the  WDF;  these  examples  can  be  evaluated  analytically  in  closed  form, 
allowing  for  close  investigation  of  the  behavior  of  the  WDF,  and  as  control 
cases  on  any  computer-written  program  for  numerical  evaluation  of  the  WDF. 
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BASIC  PROPERTIES  OF  THE  WDF 

DEFINITIONS 

A  natural  definition  of  the  time-varying  correlation  of  a  nonstationary 
complex  stochastic  process  s(t)  Is 

R(t,T)  =  s(t  +  |)  s*(  t  -  f)  .  (1) 

where  the  overbar  denotes  an  ensemble  average.  The  "center*  time  In  (1)  Is 
t,  while  the  "separation"  time  Is  T.  However,  If  an  ensemble  is  not 
available,  or  if  s(t)  is  a  deterministic  waveform,  the  obvious  extension  of 
(1)  is  simply 

R(t,T)  =  s(t  +  f)  s*( t  -  f)  .  (2) 

This  quantity  is  interpreted  as  the  instantaneous  correlation  of  waveform 
s(t)  at  time  t,  for  separation  (or  lag)  T. 

The  associated  "spectrum"  at  time  t  Is  then  available,  as  usual,  by 
Fourier  transforming  (2)  on  separation  variable  T,  to  get  at  frequency  f, 

W(t,f)  =  |  df  exp(-12wfr)  R(t,T)  = 

-  j  dr  exp(-l2*fr)  s ( t  +  |)  s*(t  -  |)  .  (3) 
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(Integrals  without  limits  are  over  the  range  of  nonzero  Integrand.  Also,  it 
Is  presumed  that  s(t)  and  Its  derivatives  decay  fast  enough  to  zero  at  t  =  too 
for  all  the  Integrals  to  converge.)  This  time-frequency  function  W(t,f)  is 
called  the  Wlgner  Distribution  Function  (WDF).  It  Is  a  real  function,  even 
when  s(t)  Is  complex,  since 

W*(t,f)  =  |dT  exp(12*fT)  s*(t  +£)  s ( t  -  \)  - 


j'  du  exp(-12*fu)  s*( t  -  jf)  s ( t  +■  jf)  =  W(t,f) 


However,  It  Is  not  necessarily  positive,  as  the  simple  example  of  a 
rectangularly  gated  pulse  quickly  shows:  for 


a  for  |t|  <  T/2 


0  otherwise 


W(t,f)  =  2E 


for  |t|  <  -  .  all  f  , 


and  zero  otherwise,  where  E  Is  the  waveform  energy: 


E  =  J  dt  |s(t)|  2  =  Jal2  T  . 


An  even  simpler  example  Is  furnished  by  waveforms  with  odd  symmetry, 
s ( -t )  =  -s(t).  Substitution  In  (3)  Immediately  yields  W(0,0)  = 

-  J tfr  |s(t/2)|  2  =  -2E.  Thus  the  origin  value  of  the  WDF  is  always 
negative  for  an  odd  waveform. 
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More  generally,  when  waveform  s(t)  is  expressed  in  terms  of  its  even  and 


odd  parts  according  to 


s ( t )  =  e(t)  *  o(t)  , 


then  the  origin  value  of  the  WDF  is 


W(0,0)  =  | dr  s{T/2)  s*(-T/2)  =  2  J"  dt  s(t)  s*(-t)  = 


=  2j dt[e(t)  +  o(t)j  je*(t)  -  o*(t)]  =  2Ee  -  2EQ  , 


where 


Ee  -  (dt|e(t)|2  ,  E0  -  Jdt|o(t)| 


are  the  energies  of  the  even  and  odd  parts  respectively.  For  nonzero  t,f, 

it  can  readily  be  shown  that  the  magnitude  of  the  WDF  Is  upper  bounded  by 

2E  -  2(E  +  E  ). 

'  e  o' 
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PROPERTIES  OF  WDF 


For  s(t)  real,  it  readily  follows  from  the  definition  of  the  WDF  in  (3) 


W(t,-f)  =  W(t,f )  for  s(t)  real  . 


In  this  special  case,  it  is  only  necessary  to  evaluate  W(t,f)  for  f  >  0. 


Define  the  voltage  density  spectrum  of  waveform  s(t)  as 


S(f )  =  J*  dt  exp( -i 2trf t )  s(t)  . 


Then  an  alternative  form  for  the  WOF  in  (3)  is 


,f)  =  J  dr  exp(-i2irfT)  s(t  +  f)  s*(t  -  f)  = 
=  J  dw  exp( i2irwt)  S(f  +  J)  S*(f  -  £)  , 


in  terms  of  S(f ) 


Integration  on  (10)  immediately  yields  the  marginals 


f  dt  W(t,f)  =  |S(f)|  2  , 
Jdf  W(t,f)  =  |s(t)J  2  , 


r  v.  >.  < 
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where  we  used  the  result 


J' dx  exp{ i2wxy)  =  6 ( y )  . 


The  quantity  in  (11)  is  the  energy  density  function,  while  that  in  (12)  is 
the  instantaneous  power.  If  we  complete  the  integrations  on  the  remaining 
variables  in  (11)  and  (12),  they  both  yield 


dt  df  W(t,f)  =  E  = 


J  dt  |s(t)|  2  =  Jdf  |S(f)|  2  . 


where  E  is  the  total  waveform  energy. 


If  waveform  s(t)  satisfies  a  time-limited  restriction,  namely 


s(t)  *  0  only  for  t1  <  t  <  t?  , 


then  (3)  reduces  to 


r 

;t,f)  =  J  dr  exp(-i2irfT)  s ( t  +  s*(t  -  f)  for  t]  <  t  <  t2  ,  (16) 


and  zero  otherwise,  where 


T  =  2  min(t„  -  t,  t  -  t,)  for  t.  <  t  <  t,  . 
m  2  I  1  c. 


Thus  the  WDF  is  time-limited  if  waveform  s(t)  is  time-limited;  however,  if 
there  are  gaps  in  s(t),  the  behavior  of  the  WDF  is  more  complicated,  as  will 
be  demonstrated  later. 
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PROOUCT  AND  CONVOLUTION 


If  waveform  s(t)  is  the  product  of  two  other  waveforms. 


s(t)  =  a(t)  b(t)  . 


then  the  WDF  of  s(t)  is  (inserting  subscripts  as  needed) 


W$(t,f)  =  JdT  exp(-12irfr)  Rs(t,r) 


=  Jdr  exp(-i2irfr)  Ra(t,r)  Rb(t,r)  = 


=  jdv  Wa(t.v)  Wb(t,f  -  v)  = 


=  Wa(t,f)  ®  Wb(t,f) 


which  is  a  convolution  on  frequency  f,  for  fixed  t. 


In  a  similar  fashion,  if  s(t)  is  the  convolution  in  time,  of  two  other 


waveforms , 


s(t)  =  a ( t ) 


then  the  WDF  of  s(t)  is 


©  b(t)  =  J 


dr  a(r)  b(t  -  r)  , 


I 


.  JP  k*  Ja* 
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Ws(t,f)  -  wa(t,f)  ®  Wb(t,f)  =  J  dr  wa(r,f)  wb(t  -  r.f)  , 


which  is  a  convolution  on  time  t,  for  fixed  f. 


AMBIGUITY  FUNCTION 


The  WDF  is  closely  related  to  the  complex  ambiguity  function  of  s(t), 
which  is  defined  here  as  [4;  section  7.2] 

X(V.T)  =  Jdt  exp( -i2trut)  s ( t  +  f)  s* ( t  -  \)  = 


■  I 


dt  exp( -i2*wt)  R(t,r)  = 


-  J  df  exp( 12«f T)  S(f  +  \)  S*(f  -  £) 


In  fact,  the  two  are  double  Fourier  transforms  of  each  other: 


dt  df  exp(i2«Tf  -  12*vt)  W(t,f) 


=  ^  dt  df  exp(i2*tf  -  i2*wt)  J  dT-|  exp( -i2irfr1 )  s(t  +  -^)  s  (t  - 

-  (J dt  dr]  exp(  -i2*»t)  s(t  +  s*(t  -  4(T  -  ^ )  = 

-  ^  dt  exp( -i2*wt)  s ( t  +  s*(t  -  J)  =  'X(v.r)  . 
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Here  we  used  (3),  (13),  and  (22) 


In  a  similar  fashion,  the  following  (single)  Fourier  transform 
relationships  on  the  WDF  hold: 


Jdf  exp( i 2irfr)  W(t,f)  =  s(t  +  |)  s*( t  -  j)  =  R(t,T)  , 


Jdt  exp(  -i2irut)  W(t,f)  =  S(f  +  ~)  S*(  f  -  j)  =  A( v,f )  . 


These  properties  are  summarized  in  the  following  diagram,  where  an  arrow 
denotes  a  Fourier  transform: 


R(t,T) 


W(t,f) 


%(«,r) 


A(v.f) 


Not  every  function  of  t,f  is  a  (legal)  WDF;  in  fact,  from  (24)  there 


follows 


Jdf  exp  ^i 2irf  ( t.j  -  t2))  W  ^  ,fj=  s(t1)  s*(t2)  =  R^1  *  -  ,  t,  -  t?j  . 


Thus,  in  order  for  a  candidate  function  W(t,f)  to  be  a  WDF,  the  function 
resulting  on  the  right-hand  side  of  (25)  must  be  separable  in  the  variables 


t^  and  y  When  and  only  when  that  separability  occurs,  the  waveform  s(t) 
can  be  recovered  from  correlation  R  or  W  (within  a  constant  unknown  phasor) 
as  f ol lows :  let 

s(y  =  Js(to)|  exp(ie(to))  ,  R(to,0)  =  |s(tQ)| 2  ,  (26) 


where  t  is  arbitrary,  except  that  s ( t q )  *  0.  Then,  from  the  right-hand 
side  of  (25), 


s(t)  = 


exp(ie(tQ)) 


for  all  t  .  (27) 


The  special  case  of  tQ  *  0  was  given  in  [3;  (17)].  The  fact  that  the 

constant  phase  e(tQ)  is  irretrievably  lost  in  R  and  W  can  easily  be  seen 

by  considering  s(t)  =  c  g(t),  for  which  W  (t,f)  =  JcJ2  W  (t,f). 

s  9 

The  box-like  function  rect(t/T)  rect(f/F)  =  1  for  |t|  <  T/2  and 

|f|  <  F/2,  zero  otherwise,  which  was  employed  for  smoothing  in  [1],  is  not  a 

WDF,  since  the  transform  on  the  left-hand  side  of  (25)  yields 

F  sinc(Ft1  -  Fy  for  jt^  +  t^j  <  1,  which  is  not  separable  in  and  y 

2  2  2  2  2 

Also,  the  Gaussian  function  exp( -t  /o  -  b  4*  f  )  is  a  legal  WDF  if  and  only 
if  b  =  o,  in  which  case  s(t)  =  (4 «o2)  1/4  exp(-t?/(2o2) ) ,  with  a  arbitrary. 


FIRST  MOMENTS  OF  W 


The  marginal  integrals  of  W  were  given  in  (11)  and  (12).  The 
(conditional)  first  moment  of  W,  with  respect  to  frequency,  is 

J  df  f  W(t,f)  =  Jdf  f  |dr  exp(-12irfT)  s(t  4-  ^)  s*(t  -  |)  - 
-  J  dt  s ( t  +  J)  s*( t  -  |)  J  df  f  exp( -i2irfT)  - 

«  £  JdT  s(t  +  |)  s*(t  -  f)  A((t)  -  Im{s*(t)  s*(t)]  .  (28) 

Here  we  used  the  result 

i2*Jdx  x  exp(i2irxy)  =  6 1  ( y )  ,  (29) 

obtainable  directly  from  (13)  by  taking  a  derivative  with  respect  to  y. 


Therefore  the  "frequency  center  at  time  t"  of  waveform  s(t)  is  defined  as 


vf(t) 


Jdf  f  W(t,f)  ^  j  Imfs'(t)s*(t3 
Jdf  W(t,f)  2*  |s(t)j 2 


(30) 


upon  use  of  (28)  and  (12).  If  we  let  complex  waveform  s(t)  be  represented 
in  terms  of  its  amplitude  and  phase  modulations  according  to 


s ( t )  =  M(t)  exptie(t)]  , 


(31) 
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then  (30)  yields  simply 


vf(t)  =  ^  e'(t)  , 


which  is  independent  of  amplitude  modulation  M(t).  (32)  can  also  be 

interpreted  as  the  instantaneous  frequency  at  time  t  of  waveform  s(t). 

Ihe  "time  center  at  frequency  f"  follows  in  an  analogous  fashion  as 

,  f .  _  -fdt  1  W(t,f)  =  _  _±  Imfs1  (f )  S*(f)) 

Pl  Jdt  W(t,f)  2ir  |s(f)|  2 


S(f)| 


If  we  represent  the  voltage  density  spectrum  S(f)  in  terms  of  its  magnitude 
and  phase, 

S(f)  -  A(f)  exp[-ia(f)]  ,  (34) 

then  (33)  reduces  to 

ut(f)  =  2*  •'(f)  •  (35) 

which  is  independent  of  A(f).  (35)  can  also  be  interpreted  as  the  group 

delay  at  frequency  f  of  waveform  s(t). 

The  unconditional  first  moments  of  W  are  frequency  center 


I**:, 


gw 
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and  time  center 

JJdt  df  t  W(t,f)  Jdt  t  |s(t)|  2 
JJdt  df  W(t,f)  Jdt  js(t)|  2 

(36)  follows  directly  from  (11)  and  (34),  while  (37)  follows  directly  from 
(12)  and  (31).  Thus,  f  is  independent  of  0(f),  and  t  is  independent  of  e(t). 


Jdt  t  M2{ t) 

J  dt  M2 ( t) 


(37) 


Alternative  forms  to  (36)  and  (37),  in  the  complementary  domains,  are 
avai lable: 


and 


Jdt 

Im{s'(t)  s * ( t )] 

l  ,  J 

f 

dt  1 

M2(t)  e’(t) 

J 

"dt  |s(t)|  2 

2ir 

*  p  » 

dt  Md(t) 

1  - 

df 

Im' 

S'(f)S*(f)}  1 

[df 

A2(f)  0’ ( f ) 

2*  J 

f  « 

)s(f )| 2  2ir 

J 

[ df  A2( f ) 

(38) 


(39) 


The  result  in  (38)  follows  from  the  use  of  (28)  and  (12)  in  definition  (36); 
a  similar  procedure  yields  (39).  The  frequency  center  f  in  (38)  is  an 

average  of  the  instantaneous  frequency  pf(t)  in  (32),  weighted  according 

.  2  — 
to  the  magnitude-squared  waveform,  M  (t);  similarly,  time  center  t  in  (39) 

is  a  weighted  average  of  ut(f)  in  (35). 
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StCOND  MOMtNIS  Oh  W 


By  taking  two  partial  derivatives  with  respect  to  T  in  (24),  there 


f  o 1  lows 


df  f2  W( t , f )  =  ~ 

8ir 


J/s  •  ( t >J  2  -  Rejs"(t)  sV,}]  - 


When  we  then  employ  (12)  and  (31),  the  (conditional)  second  moment  with 
respect  to  f  develops  into  the  form 


.  niti" 


[df  W(t,f)  8ir2  |_M2( t) 


Mm  .  r©itii 

M(t)  [_  2ir  J 

J 


Therefore  the  instantaneous  "mean-square  frequency  spread"  is 


of(t)  a  - 


Jdf  [f  -  pf(t)]2  W(t,f) 


[df  W(t,f) 


1  M  (t)  M(t) 
8w2  M2( t) 


J _ d 

,  2  dt 


(Mill] 

jM(t)  j  ’ 


where  we  employed  (32)  and  (12).  This  result  does  not  depend  on  phase 
modulation  e(t).  However,  it  should  be  observed  that  this  quantity  can  be 
negative;  consider  the  example  M(t)  =  exp(  -tv)  for  t  >  0,  with  0  <  v  <  1. 
Thus  (42)  can  not  be  interpreted  as  a  true  variance.  This  unfortunate 
feature  of  the  WDF  is  due  to  the  fact  that  W(t,f)  can  go  negative  for  some 
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The  unconditional  second  moment  with  respect  to  f  follows  from  (40)  and 
(11),  respectively,  as 

JJdt  df  f2  W(t,f)  =  ~  Jdt  |s‘(t)|  2  =  Jdf  f 2  |s(f )|  2  .  (4: 


Analogous  relations  for  the  second  moments  with  respect  to  t  can  also  be 
derived  via  a  similar  approach. 


M0MLN1S  OF  W 


The  marginal  integral  of  the  square  of  the  WDF  with  respect  to  f  is,  vi 
(3),  (13),  and  (2), 

Jdf  W2(t,f)  -  Jdf  JJ  dt  du  exp(-i2*f(T  -  u))  R(t,T)  R*(t,u)  = 


=Jdr|R(t,T)| 2  =  Jdr  Js(t  +  f)| 2  |s(t  -  f)| 2  = 

"  2[|s(t)|2®  |s(r)|2]  , 

L  Jr  =  U 

2 

which  is  the  convolution  of  Js|  with  itself,  at  argument  2t. 
complementary  result,  integrating  with  respect  to  t,  is 


f dt  w2(t,f)  =  fd»  |s( f  +  5)1 2  k(f  -  d 


If  we  complete  the  integrations  in  (44)  and  (45)  on  the  remaining 


variables,  both  yield  the  result 


dt  df  W2(t,f)  =  E2  ; 


(46) 


see  (14).  Also  note,  for  comparison,  that  the  double  integral  on  W  yielded 
E. 


Although  the  results  in  (44)  and  (45)  are  not  overly  simple,  continued 
integration  does  yield  a  surprisingly  simple  result;  multiplying  (44)  by  t, 
there  follows 


jjdt  df  t  W2(t.f)  =  Jdt  t  2  J  dx  |s(x)|  js(2t  -  x)| 


jdx  |s ( x ) |  2  J  dt  t  |s(2t  -  x )|  =  J  dx  |s(x)j  J  dy  |s(y)|  (y  +  x)/2  = 


-  \  Jdx  |s ( x )|  [tE  4-  xE]  =  l[t  E2  +  t  E2]  =  E2  t  . 
Here  we  used  (44),  (37),  and  (14).  Ihus 


Jjdt  df  t  W2(t,f) 

r  2 
jdt  t  |s(t)| 

JJ  dt  df  W2(t,f) 

Jdt  Js(t)l  2 

(47) 


(48) 


from  (47),  (46),  and  (37).  This  result  In  (48)  is  the  same  as  (37),  but  now 
2 

for  W  rather  than  W. 
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.2  • 


Conditional  first  moments  of  W  with  respect  to  f  and  t  are  also 
deriveable;  for  example, 


[  df  f  W2(t.f )  =  ^  [dr  R*(t,t)  R(t,r)  - 


af 


-  ^  f  dr  |s(t  +  |)|  2  lm|s'(t  -  |)  s*(t  -  |)]  - 


=  J  J*  dx  |s(2t  -  x)|  2  Imjs'(x)  s*(x)J  = 


-  ~  [dx  M2(2t  -  x)  M2(x)  © 1  ( x )  . 


(49) 


Here  we  used  (3),  (29),  (2),  and  (31).  When  normalized  by  the  quantity  in 
(44),  the  result  is  considerably  more  complicated  than  the  corresponding 
result  for  W  in  (30)  and  (32);  nevertheless,  continued  integrations  simplify 
tremendously.  In  particular,  there  follows,  from  (49),  (14),  (38),  (46), 
and  (36), 

2 

(50) 


Sf dt  df  f  W2( t , f )  Jdf  f  |s(f)l 

[[  dt  df  W2(t,f)  Jdf  |s(f)|  2 


This  is  the  dual  relation  to  (48),  but  derived  by  means  of  a  different 


approach.  Comparison  of  (50)  with  (36)  reveals  the  same  result  for  W  as  for 


«2. 
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CROSS  W1GNER  DISTRIBUTION  FUNCTION 


The  cross  WDF  of  two  complex  waveforms  a(t)  and  b(t)  is  a  generalization 


of  (3)  and  (10)  according  to 


Wab(t,f)  "  [dr  exp(-i2»fT)  Rab(t,r)  = 

=  [dT  exp(-i2*fr)  a ( t  +  *)  b*(t  -  J)  = 
=  Jd«  exp( i2ir«t)  A ( f  +  |)  B*(f  -  |)  , 


which  is  generally  complex.  If  a(t)  and  b(t)  are  nonzero  only  for 
a-!  <  t  <  a?  and  ^  <  t  <  b^,  respectively,  then  the  integral  limits 
on  X  in  (51)  are  explicitly  T, .  r?.  where 


Tj  =  2  max  (a1  -  t,  t  -  b?)  ,  =  2  min  (a?  -  t,  t  -  b1 ) 


If  >  X^,  then  Wab  is  zero. 


The  following  properties  of  the  cross  WDF  result  immediately: 


=  Wab<l’-f> 


wha(t.f)  -  W  .  (t,f)  , 


f  df  W,K(t.f)  =  a(t)  b*(t)  , 


ii'A* 
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[ dt  Wab(t,f)  =  A(f)  B*(f )  , 

[[dt  df  Wab(t,f)  =  [ dt  a ( t >  b*(t)  =  Jdf  A(f)  B*(f )  , 


[jdt  df  |wab(t.f)|2  -  Ea  Eb  , 


df  waa(t,f) 


wbb(t-f>  '  |f' 


2 

dt  a ( t )  b  (t) 


=  lff 


dt  df  wab(t.f)  , 


[jdt  df  Wab(t,f)  w*d(t,f)  =  [dt  a ( t )  c*( t) *  [dt  b*(t)  d(t)  . 


The  last  three  relations  follow  upon  substitution  of  (51),  interchanging 
integrals,  and  the  use  of  (13).  Again,  the  double  Fourier  transform  of  the 
cross  WDF  yields  the  cross  ambiguity  function: 


[[  dt  df  exp(-i2irvt  +  i2wfr)  Wflb(t,f)  = 


=  [  dt  exp(-i2irvt)  a ( t  f  |)  b  (t  -  ^)  = 


df  exp( i2wft)  A(f  +  f)  B  (f  -  |)  = 


=  ^ab(w,T)  =  [  dt  exp(-i2«vt)  Rab(t,r) 


S’, 
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lhe  magnitude  squared  cross-WDF  of  waveforms  a  and  b  is  also  related  to 
the  auto-WDFs  by  means  of  a  double  convolution: 


Wab(t*f)l  2  a  Jf^l  <*2  exp(-i2*fCT1  -  T2))  * 


* 


'  «  e«P( -i2.fr)  «aa(t  »  f  .  *)  Rbb(t  -  f  .  4)  - 

-  J|dt'  dx  exp(  )  Raa(t  *  |  .  Jdf'  exp^2.f'  Wbb  ^  -  f,  f)  ■ 

‘  JI111  df-  wbb(l  -  f  •*')  J«'  «p(-^*(2f  -  f ■>  *)  «aa(t  *f  . 


■ 2  wbb(i  -?  •')  waa(i  *!'v  - f)  ’ 

■  ft*  d-Haa(t  *!  • f  *!Hb(  -F  - f  ?)• 

where  we  let  T=  (t^  t-T,,)/2.  T'  -  in  the  third  line. 


( 

l 

(54)  j 

i 

t 
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NARROWBAND  REAL  WAVEfORM 


When  waveform  s(t)  is  narrowband  and  real,  it  can  be  expressed  in  terms 
of  its  low-pass  complex  envelope  c ( t )  according  to 


s ( t )  =  2  Re{c(t)  exp(i2«fQt)} 


=  c(t)  exp(i2irfQt)  *-  c*(t)  exp( -i2«fQt)  , 


where  f  is  the  carrier  or  center  frequency  of  s(t).  The  WDf  of  s(t)  is 
then  expressible  as 


Wss(t'f)  =  1  dT  exP('i2irft)  s(t  *  f)  sV  -  f>  = 


=  wcc<t'f  ‘  V  *  Wcc(t‘  '  f  -  fo>  * 


f  2  Refwcc*(t,f)  exp(i4*fot)} 


Here,  we  substituted  for  s(t),  and  used  (51)  and  (52).  Since  complex 
envelope  c(t)  is  low-pass,  a  representative  contour  plot  of  (55)  appears  as 
shown  in  figure  1.  The  wiggly  lobe  centered  at  f  -  0  is  subject  to  rapid 
oscillations  in  t,  whereas  those  lobes  centered  at  ±f  are  slowly  varying 
with  f  and  t.  A  small  amount  of  averaging  in  time  would  wipe  out  the 
undesired  oscillating  lobe,  but  maintain  the  desired  components  at  f  -  if  . 


Figure  1.  WDF  for  Narrowband  Real  Waveform 


SAMPLING  PROPERTIES 

By  letting  u  =  Z/2  in  (3),  the  WDF  becomes 

W(t,f)  =  2  ^ du  exp(-i4«fu)  s(t  +  u)  s*(t  -  u)  ,  (56) 

where  we  again  now  allow  general  complex  s(t).  If  this  integral  is  to  be 
evaluated  numerically  on  a  computer,  we  will  need  to  sample  the  integrand  at 
some  increment  A^ ,  and  apply  some  integration  rule.  In  particular,  if  we 
use  the  Trapezoidal  rule  and  carry  out  the  summation  over  -«o,  +-<*>,  we  have 
approximation 

W(t,f )  s  2A^  ^  exp(  -i4«f kA^)  s(t  +•  kA^)  s  ( t  -  kA^) 
k 

for  all  t.f.  Since  it  is  immediately  seen  from  (57)  that 


A.  -1 

it  follows  that  W(t,f)  has  period  (2At)  in  f,  when  waveform  s(t)  is 
sampled  at  increment  A^.  In  fact,  it  can  be  shown  that  W  is  the  aliased 
version  of  W: 

W(t,f)  =  ^>  W(t,f  -  ~)  .  (59) 

n 

Thus,  W(t,f)  need  be  evaluated  only  over  one  period,  say  (0,.5/a  ). 

Since  (57)  cannot  be  evaluated  for  all  continuous  values  of  t  and  f, 
we  will  limit  its  evaluation  to 


t  -  mat  ,  f  2N rA 


f“t 


where  m,  n,  Nf  are  integers.  Then  (57)  becomes  (exactly) 


(60) 


w(mAt,  2^—- j  =  2At^exp(  -i2«nk/Nf )  s((m  *■  k)A^  s*((m  -  k)At)  ,  (61) 

f  1  k 

the  right  side  of  which  is  recognized  as  an  Nf-point  discrete  Fourier 
transform.  If  the  number  of  nonzero  samples  in  k  is  greater  than  Nf,  we 
simply  collapse  them  mod  Nf,  without  loss  of  accuracy,  see  [5;  page  7]. 

A#  —  ] 

Since  the  period  of  W(t,f)  is  (2A^)  in  f,  we  only  need  consider 
0  <  n  <  Nf-1,  that  is,  0  <  f  <  (2A^)  \  Values  of  m  must  be  considered 
wherever  the  summand  of  (61)  is  nonzero. 


A  plot  of  two  of  the  infinite  number  of  lobes  of  W(t,f)  in  (59)  is 
depicted  in  figure  2  for  a  representative  bandpass  analytic  waveform  s(t). 
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The  spreads  of  the  desired  WOF  term  W(t,f)  are  7  and  B  in  time  and  frequency, 
respectively.  In  order  to  guarantee  that  aliasing  is  insignificant  in 
figure  2,  we  must  choose 


( 2At ) >  B,  that  is,  At  <  (2B) 


-1 


(62) 


1 


For  N^.  equal  to  a  power  of  2  in  (61),  an  FF1  can  be  employed  to  evaluate  W 
and  will  give  the  vertical  slice  in  f  indicated  in  figure  2  between  f  -  0  and 
f  -  ( 2Aj. )  1 ,  for  the  particular  m  value  under  consideration.  Since  the 
spacing  of  frequency  values  in  (61)  is  (2NfAt)  1 ,  then  in  order  to 
keep  track  of  the  wiggles  in  W(t,f)  as  a  function  of  f,  we  must  choose 


(NfAt)_1  <  T'1,  that  is.  Nf  >  T/At  >  2B1 . 


(63) 


Thus  the  FFT  size  may  have  to  be  quite  large  for  an  extended  WDF  in  t,f 
space. 

If  s(t)  is  real,  then  (8)  applies,  meaning  that  W  in  (61)  need  only  be 
computed  for 


0  <  n  <  N^/2,  that  is,  0  <  f  <  (4a^)  1  . 


(64) 


The  pertinent  approximate  WOF  W  is  depicted  in  figure  3.  In  order  to  avoid 
aliasing  now,  we  must  have 


(4At)  1  >  fH,  that  is,  At  <  (4f H)  1  , 


(65) 


1, 


26 
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where  fu  is  the  highest  frequency  contained  in  s(t).  This  sampling  rate 
is  twice  as  fast  as  the  usual  Nyquist  rate  for  waveform  s(t),  and  is  due  to 
the  unavoidable  factors  of  1/2  in  definition  (3). 


The  procedure  described  above,  in  (61)  et  seq.,  realizes  a  slice  in  f, 
at  fixed  t,  of  the  WDI-;  see  figures  2  and  3.  An  alternative  procedure  for 
obtaining  slices  of  the  WDF  in  t,  at  fixed  f,  is  described  in  appendix  A; 
however,  starting  with  time  samples  of  s(t),  it  requires  an  additional 
large-size  FF1  to  start  the  calculations. 
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EXAMPLES  OF  WDF 


In  this  section,  we  present  several  examples  of  the  WDF  for  waveforms 


that  are  likely  to  be  encountered  in  practice,  and  that  are  amenable  to 


simple  closed  form  solution.  A  significant  shortcut  in  the  presentation  is 


possible  when  it  is  observed  from  (3)  that  if 


r(t)  =  s ( t  -  tQ)  exp(i2irfQt  *-  i©Q)  , 


which  corresponds  to  a  time  delay  and  frequency  shift,  then 


W  (t.f)  r  h  (t  -  t  ,  f  -  f  )  . 
r  s  o  o 


Thus  we  can  choose  any  convenient  origin  for  the  waveform  s  in  time  and 


frequency,  without  loss  of  generality,  and  then  merely  shift  the  WDF 


according  to  (67),  as  appropriate. 


We  will  place  heavy  emphasis  here  on  combinations  of  Gaussian  pulses. 


both  because  of  their  analytic  tractability  and  due  to  the  fact  that  any 


waveform  can  be  expanded  into  elementary  waveforms  consisting  of  Gaussian 


wavelets;  see,  for  example,  Gabor's  original  paper  [6,  part  1,  section  5]. 


In  the  following,  frequent  use  will  be  made  of  the  following  integral: 


J  dx  exp(-  1  ax2  ±  Bx)  -  exp(lb) 


for  c*r  >  0  , 
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where  a  and  (3  can  be  complex,  with  components 


a  -  a  C  ia,  ,  (3  -  (3  +•  i  (3  .  . 

r  1  r  i 


Also,  as  a  special  case,  there  follows 


fdx  exp( -  ^  ax2  t  Ox ) 


.  a  (02  -  (32)  +  2a. [3  0. 

2ir  rv  r  l  i  r  i  /7nx 

,  1/2  P  2,2  ’  (?0) 

(a2  4-  tt2)  L  ar  +  ai  J 

'  r  i 


written  out  in  terms  of  purely  real  quantities. 


GAUSSIAN  WAVEFORM 


Let  waveform 


ao  exp(-  J  ;  a0  complex 


(Parameters  will  be  real  unless  indicated  otherwise.)  Use  of  (3)  and 


(68)  yields  WDF 


W(t,f)  -  2E  exp  - 


[-  4  -  o-f-/  • 

L  °o  J 


where  E  is  the  waveform  energy: 


E  -  v? h„l 2  •>„ 


lhe  WDF  consists  of  a  single  positive  lobe  in  t,f  space,  centered  at  the 


origin . 


■v 


kA 


S 


■a 

& 
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Observe  that  W(0,0)  is  equal  to  2E  for  this  example;  in  fact,  from  (3), 


W(0,0)  -  J  dr  s(r/2)  s*(  -T/2)  =  2E  if  s(-t)  =-  s(t)  . 


Thus  waveforms  s(t)  with  this  even  symmetry  result  in  peak  WDF  values  of  2E 
at  the  origin.  However,  if  s  has  odd  symmetry  about  zero,  s(-t)  =  -s(t), 
then  W(0,0)  =  -2E. 


The  contours  of  equal  height  of  the  WDF  in  (72)  are  ellipses.  The 
contour  for  the  case  where  the  levels  are  down  to  exp(-l)  of  their  peak 
value  is  the  ellipse  indicated  in  figure  4.  The  area  of  this  particular 
level  ellipse  is  1/2  in  the  t,f  plane.  When  this  area  is  multiplied  by  the 
peak  height  of  2E,  the  product  is  E,  which  is  just  the  volume  under  the  WDF; 
see  (14).  Thus  the  "effective  extent"  of  the  WDF  in  (72)  is  that  given  in 
figure  4,  for  relative  level  1/e  of  the  peak. 


Figure  4.  Contour  of  WDF  (72)  at  1/e  Relative  Level 
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GAUSSIAN  -MODULATED  10NE 


s(t)  =  bQ  cos ( 2irf Qt  *  0Q)  exp  J. 


The  energy  of  this  waveform  is 


E  =  2^bo  °o[1  f  cos(2*o)exp(-y2)]  * 


and  its  WOE  is 


^bo  °0)  w(t’f)  =  exp[-x2  -(y  -  yQ)2]  +  expfx2  -(y  +■  yQ)2] 


+  2  cos(4irfQt  +■  20q)  exp[-x2  -  y2]  , 


where  dimensionless  variables 


x  =  t/aQ  ,  y  =  2irf oQ  ,  yQ  =  2vfQcQ  . 


There  are  two  positive  lobes  centered  at  (t,f)  =  (0,fQ)  and  (0,-fQ),  each  of 
peak  height  approximately  E  (if  yQ  is  large).  The  contours  of  each  of  these 
lobes  are  circles  in  the  x,y  plane,  or  ellipses  in  the  t,f  plane,  as  indicated 
in  figure  4.  There  is  also  an  oscillating  lobe  centered  at  the  origin;  this 
is  an  example  of  the  general  situation  depicted  in  figure  1. 


It  should  also  be  observed  from  (77)  that  if  a  slice  in  frequency  is  taken 
of  the  WDF,  at  fixed  time  t,  that  there  is  no  fast  oscillation  in  any  of  the 
three  lobes.  Whatever  value  of  the  cos  is  encountered,  that  value  is  main¬ 
tained,  and  the  only  variation  with  y  (frequency)  is  the  Gaussian  dependence. 
Thus  if  we  locally  averaged  the  WDF  with  respect  to  frequency  alone,  that 
would  not  eliminate  the  undesired  oscillating  lobe  centered  at  (0,0). 


I 
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MULTIPLE  MOOULATED  TONES 


Consider  complex  waveform 


s(t)  = 


f  (t  -  \)2~]  ,  , 

^>ak  exp  i2irf kt - => —  ;  jakJ  complex  . 

k  L  2ok  J 


This  is  a  collection  of  tone  bursts  centered  at  ftk»fk}  in  the  t,f 

2 

plane,  with  energy  Ek  =  Yfl*Jak|  °|<-  The  corresponding  WDf  follows 


from  (3)  and  (68)  as 


W(t.f)  =  22  Ek  exp 
k 


+  4f?  a.  exp 


L  4 

~  (t  -  s 
_2 

L  °tt 


~  -  J  4ir2(  f  -  f.  )2  * 


-  ^  4w2(f  -  f 


.>•] 


r  *  r  ~2  a  -  h  1  -  h 

Rejaka^  exp 2ir( f R  -  fg  )t  +  2*(f  -  f^)^- — 1 - - 2~ 


where 


\  =  26k  *  ti)  •  fm  ‘  ?(fk +  !*)  - 

-2  l/2  .  2\  J_  _  l/l  .  _l\ 

^  x)’  V  °2J 


i; 


Ihe  first  line  of  (80)  represents  the  desired  positive  lobes  centered  at 


(t^.f^),  each  scaled  according  to  its  energy.  The  remaining  undesired 
lobes  are  centered  at 


\  *  V 


Vi 

2 


for  all  k  i t  Jl  , 


(82) 


and  oscillate  with  t  and/or  f.  lhese  locations  in  (82)  are  halfway  between 
every  possible  pair  of  desired  lobes;  their  amplitudes  are  proportional  to 
the  geometric  means  of  the  corresponding  interacting  lobes,  and  therefore 
constitute  significant  interference  effects  to  interpretation  of  the 
computed  WDF.  Furthermore,  the  locations  in  (82)  can  occur  in  time  where 
the  waveform  s(t)  is  zero,  and/or  in  frequency  where  the  spectrum  S(f)  is 
zero.  This  most  undesirable  feature  of  the  WDF  has  been  reported  previously 
in  [7,8].  The  only  saving  feature,  that  should  allow  salvaging  the  WDF,  is 
that  the  undesired  lobes,  k  < Jl  in  (80),  oscillate  posit ive-and -negative  and 
can  be  averaged  out  by  smoothing  the  WDF.  Of  course,  via  this  smoothing 
procedure,  the  desired  lobes  will  also  be  smeared  somewhat,  but  this 
trade-off  appears  to  be  required  in  order  to  make  a  meaningful,  useful 
interpretation  of  the  WDF  at  all  points  of  the  t,f  plane. 

Ihe  envelope  of  the  k,Jl  lobe  in  (80)  is  proportional  to  an  exponential 


of  an  elliptical  function.  When  this  exponential  has  decreased  to  1/e  of 
its  peak  value,  the  corresponding  elliptical  contour  has  area 


in  the  t,f  plane,  the  latter  value  of  1/2  being  the  area  of  every  desired 
lobe.  Thus  the  undesired  oscillating  lobes  are  smeared  out  more  than  the 
desired  positive  lobes. 

If  we  restrict  (79)  to  two  equal-duration  bursts  with  the  same  time 
center,  but  different  center  frequencies,  the  undesired  lobe  oscillates  only 
with  t,  not  f.  This  is  similar  to  example  (75)-(77).  On  the  other  hand,  if 
(7^)  is  restricted  to  two  equal -duration  bursts  with  different  time  centers, 
but  the  same  frequency  center,  the  undesired  lobe  oscillates  only  with  f, 
not  t. 

More  generally,  for  two  equal -duration  bursts  with  different  time  and 
frequency  centers,  the  undesired  lobe  has  no  fast  oscillation  along  lines  in 
the  t,f  plane  which  are  parallel  to  the  line  joining  the  centers  of  the  two 
positive  lobes  in  the  WDF .  For  two  unequal -duration  bursts,  the  situation 
is  more  complicated,  and  there  is  generally  oscillation  along  all  straight 
lines  in  the  t,f  plane. 

What  these  simple  examples  demonstrate  is  that  if  we  want  to  locally 
smooth  (average)  the  WDF,  in  an  effort  to  wipe  out  the  undesired  oscillating 
cross -terms,  that  smoothing  must  be  applied  in  both  t  and  f,  not  either  one 
alone.  Of  course,  such  smoothing  will  also  tend  to  smear  the  desired 
positive  lobes;  thus  the  minimum  amount  of  smoothing  to  guarantee  a 
nonnegative  WDF  is  of  interest. 
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Although  these  conclusions  have  been  drawn  from  the  particular  example 
of  Gaussian  modulated  tone  bursts  in  (79)  (for  analytic  simplicity),  they 
hold  generally.  Appendix  8  demonstrates  the  oscillating  character  of  the 
interacting  cross -terms  of  the  WDF  for  a  waveform  with  two  separated  energy 
bursts  in  time  of  general  shape. 


Ihe  ambiguity  function  of  waveform  s(t)  in  (75)  is  considered  in 
appendix  C.  It  has  some  similar  properties  to  the  WDF  and  some  significant 
differences,  which  make  it  much  less  desirable  as  a  descriptor  of  a  signal's 
concentration  in  time-frequency  space. 

I.1NFAR  FRF.QUFNCY  MODULATION 

Fiere,  we  consider  waveform 


s ( t )  =  aQ  exp 


t2  .  °o  .2 
~2  +  1  1  1 
-  2°o 


a  >  0  ,  a  complex  .  (84) 

o  -  o 


lhe  instantaneous  frequency,  according  to  (31)  and  (32),  is  a  linear 
function  of  time, 


wf(t)  =  2 it  1  • 


(85) 


while  the  envelope  is  Gaussian.  When 


t  =  t0  =  t  o0 Y*/2'  , 


the  magnitude  of  the  waveform  s(t)  is 


(86) 


|s(tQ)|  -  |aQJ  exp(  -w/4)  -  .456  |s(0)|  . 


(87) 

35 
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If  we  define  the  duration.  At,  of  s(t)  as  the  time  between  these 


function  values,  then 


At  -  y2«'  . 


During  this  time  interval,  the  instantaneous  frequency  in  (85)  sweeps 


through  a  bandwidth 


Af  "  2«  At  1  VoA57’- 


Ihe  t ime -bandwidth  product  of  waveform  ${t)  is  therefore 


At  Af  =  a  o  S  0  (>  0)  , 

0  0  0 


when  the  time  duration  is  defined  as  the  interval  between  the  function 
values  in  (87).  Ihis  quantity,  0q,  is  an  important  parameter  of  the 
linear  frequency  modulation  waveform  (84). 


Ihe  WDF  of  (84)  follows,  upon  use  of  (3)  and  (68),  as 


W(t,f)  =  2E  expl- 


t£ 

2 

_  °o 


a  (2wf  -  a 
0 


-  2E  exp[-x2  -(y  -x0Q)2]  = 

-  2E  exp£-x2(  1  ♦  e2)  f  2xy  0q  -  y2J  , 


where  we  employed  (78)  and  (90).  This  is  an  everywhere  positive  lobe 
centered  at  the  origin  of  time-frequency  space,  with  contours  that  are 
tilted  ellipses.  The  peak  value,  2E,  is  independent  of  the  amount  of 
frequency  modulation. 


«  V  ’  *  V 


Tor  a  given  value  of  time  t,  the  frequency  f  that  maximizes  the  WDf  in 
(91)  is 


f  2^  1  .  that  is.  y  =  xeQ  . 


(92) 


which  is  just  the  instantaneous  frequency  in  (86).  However,  this  line, 
(92),  in  the  t.f  plane  is  not  the  major  axis  of  the  elliptical  contours  of 
the  WDf.  A  similar  observation  regarding  the  ambiguity  f unct ion  ,T) , 
(23),  of  the  linear  frequency  modulation  waveform,  namely 


X(Y,T)  =  t  exp[-  j{x>2(l  .  e2)  -  2 x ' y ' 0Q  *  y,2}j  . 


(93) 


where 


a  *  y'  *  2*v°o  * 
0 


(94) 


has  been  made  in  [4;  page  124], 


What  this  means  is  that,  if  the  WDF  of  a  waveform  is  evaluated 
numerically  from  a  given  data  sequence  (via  (61)  for  example),  then  the  tilt 
of  the  major  axis  of  the  contours  of  the  computed  WDF  is  not  directly  the 
amount  of  linear  frequency  modulation  in  the  waveform.  Rather,  the  major 
axis  of  the  ellipse  in  (91)  lies  along  the  line 


y  =  x 


tan  ir 


(95) 


in  the  x,y  plane,  where 


tan  3|z 


«■  e 

VO  0 


>  e. 


(96) 


(See  appendix  0  for  detailed  derivations  on  the  rotation  of  coordinate 
axes.)  Thus,  the  major  axis  (95)  of  the  ellipse  is  more  tilted  than  the 
instantaneous  frequency  line  (92).  Equation  (96)  can  be  inverted  and  solved 
for  the  linear  frequency  modulation  parameter  0q  according  to 

e0  -  tan^  -  1/tan-^r  .  (97) 

in  terms  of  the  measured  or  calculated  major  axis  tilt,  tan^,  in  the  x.y 
plane.  The  detailed  procedure  for  solving  for  both  oq  and  ao,  in  terms 
of  a  computed  WDf  in  the  t , f  plane,  is  discussed  in  the  example  in  appendix 
D,  especially  (D-28)  and  (D-29). 

When  the  exponential  in  (91)  is  down  to  1 /e  of  its  peak  value,  the 
ellipse  at  that  level  has  area  *  in  the  x.y  plane.  This  may  be  seen  by  use 
of  (01)  and  (0-20),  with  A  -  1  02,  B  =  -2©0,  C  =  1 ,  D  -  E  =  0,  F  =  -1, 

for  which  G  -  1  via  (0  19).  This  corresponds  to  area  1/2  in  the  t,f  plane, 
as  seen  by  (78).  Therefore,  the  peak  height,  2E,  times  the  "effective"  area 
is  again  E,  as  it  was  for  the  simple  Gaussian  pulse  of  (71)  and  figure  A. 
Thus,  although  the  volume  of  the  W0F  in  (91)  has  been  redistributed  in  the 
t,f  plane,  by  virtue  of  linear  frequency  modulation,  the  effective  area  is 
maintained,  although  now  located  as  a  tilted  ellipse. 

A  plot  of  the  ellipse  of  (91)  at  the  1/e  level,  namely 

X2  *-  (y  ~  xe0)2  =  x2(l  f  e2)  -  2xy@0  +  y2  -  1  ,  (98) 

is  given  in  figure  5,  when  ©q  -  1.5.  The  instantaneous  frequency  line 
(92)  as  well  as  the  major  axis  (95)  are  delineated,  and  are  clearly  seen  not 
to  overlie  each  other. 
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GATED  LINEAR  FREQUENCY  MODULATION 


All  the  previous  examples  in  this  section  had  Gaussian  envelopes.  We 
now  consider  a  rectangularly  gated  waveform  with  linear  frequency  modulation: 

s ( t )  =  aQ  exp^i  t^j  for  |t }  <  |  ;  aQ  complex  .  (99; 

Equation  (3)  yields  directly  WOE 


s ini ( 2wf  -  a  t)(T  -  2  It  I ) J 


W(t.f)  -  2E  ( 2trf  -  aQt)T 


for  |t|  <  5  .  all  f  ,  (100 


and  zero  otherwise.  Along  the  instantaneous  frequency  line.  (85),  the  WDE 
is  2E  (1  -  2  \t l  /T )  for  It l  <  T/2,  which  Is  nonosci I latory  and  positive. 
However,  in  other  portions  of  the  t,f  plane,  (100)  does  go  negative,  due  to 
the  sin  term. 


For  a  given  value  of  t,  the  quantity  W  in  (100)  is  maximized  by  choosing 

f  =  oQt/(2n),  but,  again,  this  is  not  the  major  axis  of  the  contours  of  the 

2  2 

WDF.  In  figure  6,  these  contours  are  plotted  for  aQT  =  1  and  <*oT  -  10. 

In  fact,  the  contours  are  no  longer  ellipses,  although  they  tend  to  resemble 

2 

ellipses  near  the  origin,  when  frequency  modulation  parameter  aQl  is  large; 
see  the  bottom  figure,  where  the  instantaneous  frequency  line  and  the 
mountain  ridge  (curve  of  slowest  descent)  have  been  sketched. 


T 


The  ambiguity  function  of  waveform  (99)  is 


X(».x>  =  l - i 


sin[^(2iru  -  qqt) ( 1  -  |t|  )  ] 


^(2iiu  -  aQT)T 


for  Jrj  <  1  ,  a  1 1  v  , 


(101) 


and  zero  otherwise.  It  is  similar  to  the  WDf  in  (100),  but  is  spread  out 
more  in  the  v.t.  plane. 


SHORT-TERM  SPECTRAL  ESTIMATION 


l.t  |gB 


'  VJVL’VV’  VW  BT  S  »  V  V 
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Some  advantageous  features  of  the  WDF  have  been  brought  out  by  earlier 
examples,  such  as  the  concentrated  positive  lobes  in  the  t,f  plane  about 
locations  corresponding  to  obvious  bursts  of  energy.  However,  the  WDF  also 
goes  negative  in  surrounding  regions,  causing  difficulty  in  interpretation; 
see  figure  1,  (80),  appendix  B,  or  [7,8].  What  is  needed  is  some  form  of 
smoothing  of  the  WDF  so  as  to  eliminate  or  suppress  the  oscillating 
components;  however,  this  averaging  must  be  two-dimensional,  carried  out  in 
both  time  and  frequency,  for  the  reasons  presented  in  the  sequel  to  (83). 

We  now  present  one  method  of  smoothing  the  WDF,  which  guarantees  a 
non-negative  distribution  in  time -frequency . 


WEIGHTED  SPEC1RAL  ESTIMATE 


The  voltage  density  spectrum  S(f),  corresponding  to  waveform  s(t),  was 
defined  in  (9)  as  the  Fourier  transform  over  all  time.  In  order  to  bring 
out  properties  which  are  local  in  time,  a  weighting  must  be  applied  before 
transformation.  In  particular,  we  generalize  (9)  to 


yt.f)  =  | dt1  exp(-i2nft1 )  sfty  u  (t  -  y  = 


=  exp(-i2irft)  j"  df1  exp(i2irtfj)  S(f.j)  U  (f  -  f^)  , 


(102) 


& 

I 


I 

I 

if 


9 

2 
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K 
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where  weighting  u  will  tend  to  be  a  narrow  function  centered  about  its 
origin;  thus  the  weighting  in  (102)  will  accent  the  behavior  of  waveform  s 
in  the  neighborhood  of  time  t.  The  function  LI  is  the  Fourier  transform  of 
u.  The  short-term  power  spectral  estimate  (at  time  t  and  frequency  f)  of 
waveform  s,  relative  to  weighting  u,  is  then  defined  as 


SuU.n 


(103) 


See  also  [2,  p.  768] . 


The  following  symmetry  properties  of  definition  (102)  follow: 


Su(t.f)  =  Us(t,f)  exp(  -i2irf  t)  , 


|su(t’f»|2  ■  K't-f,l*  • 


(104) 


where  U$  is  the  spectrum  of  waveform  u  relative  to  weighting  s.  That  is, 
Us  is  the  dual  of  $u .  Also,  by  use  of  (53),  we  can  express 


su(t.f)  =  ‘Xs0(f ,t)  exp(-i.ft)  , 


in  terms  of  the  complex  cross -ambiguity  function  of  s  and  IT,  where  u  is  the 
mirror  image  of  u:  u(t)  =  u(-t).  Also,  the  same  shifting  property,  given 
in  (66)  and  (67)  for  the  WDF,  holds  as  well  for  quantity  (103). 


RELATION  TO  WOFs 


There  is  a  very  important  relation  between  the  short-term  spectral 
estimate  (STSE),  (103),  and  the  WOFs  of  s  and  u;  namely,  by  use  of  (102), 
( 2) ,  and  ( 3) ,  we  have 


\jn\x  fyowi  vwr/vv  v  v  w^’WAWTV  vvv  /  w.*r*  w*JW^:wwiF.nf^ 
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|Su(t,f)l  2  =  J/dt!  dt2  exp^-12*f(t1  -  t2))  s(t1)s*(t2)u*(t  -  t,  )u(t  -  t?)  - 
=  jjdT  dt1  exp(  - i 2-n-f XT )  s^t’  +  |)  s*^f  -  |)  u  (t  -  t'  +  u*  ^t  -  t1  -  |)  = 
-  JJdTdt'  exp(-i2irfx)  R^t'.x)  Ry(t  -t'.x)  = 

=  JJ  dt'  df 1  Ws(t',f)  Wu(t  -  t\f  -  f ' )  = 


=  W$(t.f)  ®  Wu(t,f)  . 


(106) 


This  relation  states  that  the  STSE  is  a  double  convolution,  in  both  t  and  f, 
of  the  WDFs  of  waveform  s  and  weighting  u.  That  is,  the  STSE  |su(t,f)J  2 
of  waveform  s,  relative  to  weighting  u,  is  a  smoothed  version  of  the  WDF  of 
waveform  s,  where  the  smoothing  function  is  the  WDF  of  weighting  u. 
Furthermore,  since  the  left-hand  side  of  (106)  can  never  be  negative,  and 
since  s  and  u  are  arbitrary,  (106)  shows  that  the  double  convolution  of  any 

two  WOFs  is  never  negative  for  any  values  of  t,f.  This  furnishes  a 

possibility  of  accomplishing  smoothing  of  a  computed  WDF  of  waveform  s,  with 

a  guarantee  of  a  nonnegative  distribution  resulting;  of  course,  W^  must  be 

a  legal  WDF,  as  discussed  in  (26)  et  seq.,  in  order  to  guarantee  this 
nonnegative  property. 


Since  Jsu|2  is  a  double  convolution  of  WDFs  W$  and  Wy,  it 
follows  that  the  double  Fourier  transform  of  the  STSE  is  given  by 


1  ,w, 
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JJdt  df  exp(  i  2irfX  -  12irtv)  ]Su(l,f)|  2  =  %u  (  » ,T)  Xs  { w  ’T)  ’  (107) 

where  ^  and  ^  are  the  complex  ambiguity  functions  of  u  and  s, 
respectively;  see  (23).  This  leads  to  an  alternative  expression  for  the 
STSE  as 

|su(t,f)|2  -  JJdt  du  e  x  p  (  —  i  2irfT  +-  i2*tv)  ^(v.r)  X.(  *.D  •  (108) 

Therefore,  if  the  complex  ambiguity  function  of  s  is  computed,  it  can  be 
multiplied  by  the  ambiguity  function  of  an  arbitrary  weighting  function  u, 
and  followed  by  a  two-dimensional  Fourier  transform.  There  is  no  need  to 
calculate  the  WDF  via  this  route;  also  several  different  weighting 
functions  could  be  utilized,  each  at  the  expense  of  a  two-dimensional 
Fourier  transform.  The  end  result  for  the  STSE  is  always  nonnegative.  Of 
course,  the  same  result  is  obtainable  directly  by  taking  the  magnitude- 
square  of  definition  (102). 

MARGINALS  OF  SPECTRAL  ESTIMATE 

There  follows,  from  (106)  and  (12),  the  marginal  relation 
J  df  |su(t,f)|  2  =  Jdt1  |s(  t  •  )|  2  |u(t  -  t*)|  2  =  |s(t)|  2  ®  |u(t)|  2  .  (109) 

2  2 

Thus,  the  time  marginal  of  |SU|*  is  not  directly  |s(t)J  ,  but  is 
smeared  by  the  weighting,  according  to  |u(t)|  2.  In  a  similar  fashion, 
from  (106)  and  (11),  the  frequency  marginal  is 


46 
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[dt  |su(t.f)| 2  =  Js( f )|  2  $  |u(f)|  2 


(110) 


Again,  Js ( f )|  2  is  smeared  by  window  |u( f )[  2. 


Mnally,  completing  either  of  the  integrations  in  (109)  or  (110),  over 


the  remaining  variable,  yields 


[fit  Of  |su(t.f)l2  =  Es  Eu  , 


(HD 


where  E$  and  Eu  are  the  energies  of  s  and  u,  respectively;  see  (14). 

Since  weighting  u  is  arbitrary  and  under  our  control,  we  can  easily  choose 
E^  to  be  1 ,  without  loss  of  generality;  then  the  volume  under  the  STSE 
will  be  equal  to  the  energy  in  waveform  s  being  analyzed,  just  as  for  the 
WOE  in  (14). 


MOMENTS  OE  SPEC1RAL  ESTIMATE 

If  we  use  (110)  and  (14),  we  find  the  following  development: 

jj  dt  df  f  |su(t,f  )\  2  =  {df  f  [  dv  |s(v)|  2  |u(f  -  »)|  2 
=  [  dv  |S(v)|  2  J  df  (f  -  v  4-  v)  |u(f  -  v)|  2  = 
-jo*  |S(*)|2[jof,  f,  [u(f,)|2  *  »EU]-- 

*  Es  j  f,K>|2  *  Eu  J"*  •  M  2  • 


(112) 


■  Vi'.'*  ,  !  .  » 


■!  r * ;•  v*  ■  *T.  •«;» 


TR  8225 


Combined  with  (111),  there  results 


Ifdt  df  f 

su(t.f)|2 

Jdf  f  Ju(f  )\ 2 

Jdf  f  |s(f)|2 

ft  «  Of 

J df  |u(f  )| 2 

J  |s(f  )|  2 

That  is,  the  first  moment  in  f  of  the  S7SE  is  the  sum  of  the  frequency 
centers  of  JllJ  and  )s|  .  This  should  be  compared  with  the 
corresponding  result  in  (36)  for  the  WOT,  where  only  the  last  term  in  (113) 
is  present.  The  presence  of  weighting  u  in  definition  (102)  adds  an 
additional  term  to  the  frequency  center  unless  |ll(f)|  is  even  about  f  -  0; 
in  this  latter  case,  (113)  reduces  to  (36). 

In  a  similar  fashion,  the  first  moment  in  t  of  the  STSt  is  found  to  be 

JJdt  df  t|su(t,f)|  2  {dt  t|u(t)l  2  Jdt  t  [s(t)|  2 
jjdt  df  |su(t.f)|2  Jdt  )u(t)|  2  Jdtjs(t)|2 

Again,  a  sum  of  time  centers  results;  but  if  weighting  |u(t)|  is  even  about 
t  =  0,  then  (114)  reduces  to  the  same  result,  (37),  as  for  the  WOF . 


i 
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CONI)  1  I  IONAL  MOMf.NI 


Just  as  in  (28) -(35)  for  the  WDF,  conditional  moments  of  the  ST SE  can  be 
defined.  For  example,  directly  from  (105),  we  have 


Jdf  f|su(t,f)|2  -  Jdf  f  |xsB(M>|2  • 


(115) 


in  terms  of  the  cross-ambiguity  function  of  s  and  u,  where  u(t)  =  u(-t). 


An  alternative  time-domain  expression  is  possible  for  the  frequency 


moment  in  (115) :  define 


g( t , t t )  «  s ( t-| )  u  (t  -  ti) 


(116) 


Then  from  (102)  and  (29), 


Jdf  f  |su(t,f)|  2  «  Jdf  f  JJ  dt1  dt2  exp  (-i2wf(t1  -  t2))  g(t,t1 )  g*(t,t2)  = 


dt1  dt2  g(t,t] )  g  (t,t?)  g  6,U1  -  t?)  = 


‘  T, 


g’  (t.t1 )  g  (t,^)  , 


(117) 


where 


g,(t’V  5  3t^  9(t’V  =  ~  {'s(t1 )  Mt  -  t^}  .  (118) 


If  we  represent  waveform  s  in  terms  of  its  magnitude  and  phase  according  to 
(31),  and  do  likewise  for  weighting  u  as 


v 

H 

•* 


u(t)  =  E(t)  exptip(t)]  . 


(119) 


then  substitution  in  (117)  results  in  the  simplified  form 

Jdf  f|Su(t.f)|2  =  ^  dt-,  I®1  ( t1 )  -  P‘(t  -  t1  )J  M2(t1)  t2(t  -  .  (120) 


When  this  result  is  combined  with  (109),  the  normalized  conditional 
first  moment  is 


[df  f|su(t,f)|2  1  Jdt]  [»'  (t] )  ~  P'(t  -  t1j]M2(t1)  E2(t  -  t^ 

Jdf  |su(t.f)|  2  2*  J  dt1  H2(t1 )  E2(t  -  t,) 


(121) 


(This  reduces  to  (38)  when  E(t)  =  1,  P(t)  =  0,  that  is,  u(t)  -  1,  in  which 
case  S  ( t ,  f )  =  S(  f ) . )  Generally,  (121)  is  an  average  of  ©'(tj)  -  P '  ( t  t] ) , 
weighted  according  to  the  instantaneous  powers  of  s  and  u. 
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EXAMPLES  OF  SHORT-TERM  SPECTRAL  ESTIMATION 


Here  we  will  reconsider  many  of  the  examples  presented  earlier  for  the 
WDF  investigation.  The  particular  example  of  weighting  u  adopted  here  in 
spectral  definition  (102)  will  be,  for  the  time  being, 


u(t) 


(122) 


where  duration  measure  a  is  under  our  control.  The  energy  E  of  this 
waveform  is  unity,  in  keeping  with  the  discussion  in  (111)  et  seq.,  which 
guarantees  that  the  volume  under  the  S1SE  will  be  the  energy  E^  -  E  of  the 
waveform  s  being  analyzed. 


GAUSSIAN  WAVEFORM 


Ihe  waveform  s  was  given  in  (71);  its  transform  Su(t,f)  is  obtained  by 


substituting  (71)  and  (122)  in  (102)  and  using  (68): 


S(t,f)  =  E  — 

U  l  a 


1/2 


a 


exp  |-  — -  if  o*  -  iirft  —  +  i  arg(ao 
j  a 

a  a 


where  E  is  the  energy  of  s  and 


(123) 


(124) 


2  2  2. 

The  quantity  a^  is  the  arithmetic  average  of  a  and  aQ,  while 

2 

is  the  harmonic  average.  The  STSt  then  follows  immediately  (or  by 
use  of  (70)  directly)  as 

1  /t?  j  2f 2 
exp  -  2\2  +  f  ° 

-  Va 

The  volume  under  ST S£  (125)  is  readily  verified  to  be  t,  as  it  must  be. 

The  half  widths  of  the  ellipse  at  the  1/e  relative  level  are  a  ,  (2*a.  ) 

3  n 

respectively,  in  the  t,f  plane.  The  area  of  this  ellipse  in  the  t,f  plane  is 


This  area  is  at  least  twice  as  great  as  that  for  the  corresponding  WDF  in 

figure  4,  and  even  then,  only  when  the  proper  guess  is  used  for  the 

weighting  u,  namely  a  -  a  .  Since  waveform  duration  a  will  likely  be 

o  o 

unknown  in  practice,  the  mismatch  factor  in  (126)  will  smear  the 
concentration  of  the  S1SE  somewhat.  For  example,  if  a  is  off  by  a  factor  of 
2  from  oq  (either  double  or  half),  (126)  is  1.25  instead  of  its  minimum 
value  of  1 . 

Ihe  area  enlargement  factor  a  /a.  in  (126)  is  also  the  same  factor 

3  n 

by  which  the  peak  of  the  STSE  in  (125)  is  down  from  its  best  value  of  E. 

Thus,  the  STSE  has  a  decreased  peak  and  enlarged  effective  area  relative  to 

the  WDF,  the  relative  factor  being  at  least  2,  and  being  a/a  +  a  /a  in 

oo 

general.  Both  distributions  contain  volume  E,  independent  of  a. 


1R  8225 


This  example  demonstrates  the  presence  of  "window  effects"  in  the  ST  St 
t  fid  t  are  not  seen  in  the  WOF .  That  is,  whereas  the  effective  ellipse  in 
figure  4  depended  only  on  waveform  parameter  aQ,  the  ellipse  here  depends 
additionally  on  weighting  parameter  a,  in  such  a  fashion  as  to  always  smear 
the  concentration  of  energy  in  the  t.f  plane  by  at  least  a  factor  of  2.  In 
trade,  we  always  have  the  guarantee  that  the  ST  St  (103)  will  be  nonnegative, 
and  that  it  will  not  contain  the  large  interference  phenomena  inherent  in 
the  WOF;  see  (80) -(81 ) . 

MULT  I PLE  MODULATED  TONES 


The  waveform  of  interest  is  given  in  (79).  The  transform  Su(t,f)  is 
found  by  use  of  (122),  (102),  and  (68): 

1/2 

exp 


SuU.f)  - 


(  -  \Y  2  2  2 

- - 2  -  * <*  -  fk>  °hk 


4  a 


ak 


i*(f  -  fk) 


2  t  .  2  t 
°k  1  *  °  tk 


ak 


+  i  arg(ak 


3' 


(127) 


where  energy  Ek  =  yw* |ak ]  «k.  and  averages 


ak 


l/2  2\  J_  l/l  A 

?  V  °k  /  ’  2  '  2  (  2  +  2  I 

°hk  ^  °k  / 


(128) 


53  I 


The  STSt  is  the  magnitude  squared  value  of  (127);  the  resulting  double 


sum  has  diagonal  terms 


°hk  1  /(*  "  lk)  ,  .  2,, 

2[krB'-;  - 2  f  4- (f 

k  ak  L  L  °ak 


( )  29) 


which  are  identical  to  (12S),  except  for  the  indexing  by  k  and  the  shift  to 
center  t^.f^  in  the  t,f  plane,  as  expected.  If  one  value  of  weighting 
parameter  a  is  used  to  evaluate  the  S1SE  for  all  t  of  interest,  it  cannot 
simultaneously  match  all  the  different  possible  values  of  for  the 

various  pulses.  This  will  cause  some  of  the  components  in  (129)  to  be  more 
severely  degraded  than  others,  in  terms  of  decreased  peaks  and  spread 
effective  areas;  the  pertinent  factor  is  again 


for  the  k  th  component  lobe.  If  some  apriori  knowledge  of  the  values  of 
is  available,  this  suggests  using  different  values  of  o  for  those 
values  of  t  near  the  corresponding  values  of  ft^} ,  in  an  effort  to 
minimize  the  factors  (130)  for  different  k. 

As  for  the  off-diagonal  terms  of  |Sy ( t , f )|  2  in  (127),  the  kJPterm  is 


proportional  to 


If  t  is  not  near  one  of  the  time  centers  {t^j  ,  or  if  f  is  not  near  one  of 
the  {fjJ.  this  term  will  be  very  small,  due  to  the  exponential  decay.  In 
particular,  halfway  between  the  dominant  desired  peaks  of  (129)  at 
the  quantity  (131)  will  be  essentially  zero.  This  is  in  distinction  to  the 
WOF  result  in  (80)  et  seq. 

LINEAR  FREQUENCY  MODULATION 


This  waveform  is  specified  in  (84).  Its  STSE  is,  upon  use  of  (70), 


1  +  r  +  r  8  )  +■ 
o 


where 


♦  y  (1  +  r)/r  -  2xy  e 


■.A- 


(132) 


8  -  a  a  ,  r 

0  0  0 


2  *  H2 
°0 


1  2 

-+2+r+re  , 

r  o  ’ 


(133) 


and  where  we  define,  here, 


x  =  -  ,  y  =  2*fo  . 
a 


(134) 


By  means  of  the  results  in  appendix  D,  the  area  of  the  contour  ellipse 
in  (132),  at  the  1/e  relative  level,  is  found  to  be 


in  the  t,f  plane. 


area 


(135) 
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Thus,  the  product  of  the  peak  height  of  the  STSE  in  (132)  and  the  effective 

area  is  again  E,  the  volume  under  the  STSE,  regardless  of  the  values  of  a, 

a  ,  e  .  (For  0=0,  (135)  reduces  to  (126),  as  it  must.) 
oo  o 


To  minimize  the  effective  area  and  to  maximize  the  peak  value  of  the 
STSE,  the  common  quantity  in  (133)  should  be  minimized.  This  is 
accomplished  by  choosing  the  weighting  parameter  a  in  (122)  as 


(1  +•  Oq) 


2  'A  ’ 

(1  +  ®q) 


(136) 


which  would  require  knowledge  of  both  the  duration  a  and  the  amount  of 

o 

frequency  modulation  8q  in  waveform  s.  Even  if  that  information  were 
available,  the  minimum  area  in  (135)  becomes 


mimnum  area  = 


1  +  1/1  +•  0 


in  the  t,f  plane  , 


(137) 


which  still  increases  as  ye^T?  for  large  8q.  Thus,  even  the  best  choice 

of  a  for  the  weighting  results  in  considerable  spreading  of  the  concentration 

ellipse  and  in  peak  reduction  of  the  STSE;  searching  in  a  is  not  overly 

helpful  because  the  simple  weighting  pulse  (122)  is  a  poor  facsimile  to  the 

linear  frequency  modulation  waveform  (84),  especially  for  large  amounts  of 

frequency  modulation,  as  measured  by  parameter  0  . 

o 


$ 

:: 

?: 

A 
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MORE  GENERAL  WEIGHTING 


lhere  is  no  need  to  restrict  the  weighting  u  in  STSE  (102)  and  (103)  to 
be  the  simple  Gaussian  pulse  in  (122).  In  this  section,  we  generalize  it  to 
allow  for  some  linear  frequency  modulation: 


u(t)  =  (*  o  )  exp  - 


[■  L*  * 1 1 1?] ; 


£u  '  '  ■ 


(138) 


The  waveform  of  interest  here  is  again  the  linear  frequency  modulation 


example 


s ( t )  -  aQ  exp  j^-  +  i  t*j  ,  aQ  >  0  , 


(139) 


as  in  (84). 


the  ST  St  follows  from  (102),  (103),  and  (70),  after  a  considerable 
amount  of  manipulations,  as 


c  ,  ,|  2  2E  f  1  f  2..  2  2  2  2 

u  ’  ■  Jp  *r.r9ttr,80)* 


4-  y  (1  +  r)/r  -  2xy(l  +  rq)e 


oil 


(140) 


where 


a  a  ..  1  „  ...22 

r=— ,  q  =  —  ,  H„=-f2+r+r(q-l  e„  , 


(141) 


i. *'•>’*.*>.**. r>,r>  .‘n1 


in  addition  to  (134).  The  quantities  r  and  q  are  mismatch  factors, 
reflecting  the  lack  of  knowledge  of  weighting  (138)  about  the  waveform  (139). 

The  area  of  the  contour  ellipse  of  (140)  at  the  1/e  relative  level  is 
(by  means  of  appendix  D) 

area  =  ^-he  t.f  plane  .  (142) 

This  is  also  the  same  factor  by  which  the  peak  of  (140)  is  down  from  E. 

Thus,  a  minimum  value  for  is  desired.  This  can  be  achieved  by  choosing 

r  -  1,  q  =  1  In  (141),  for  which  the  minimum  H3  =-  4  and  the  minimum  area  ---  1; 

however,  this  requires  that  we  choose  a  =  a  and  a  =  a  ,  which  is  not  a 

oo 

likely  situation  in  practice,  without  some  apriori  knowledge  about  the 
waveform  s.  If  this  fortuitous  situation  of  perfect  match  of  the  parameters 
does  occur,  the  STSE  in  (140)  reduces  to 

|su(t,f)|  2  =  E  exp£-  j{x2(l  +  e2)  -  2xy  Gq  +  y2j^|  ,  (143) 

which  is  identical  to  the  corresponding  WDF  in  (91),  except  for  a  factor  of 
2  outside  and  inside  the  exponential.  Thus,  the  effective  area  is  doubled 
and  the  peak  is  halved. 

As  special  cases  of  weighting  (138),  if  a  =  0  (no  frequency  modulation 
in  the  weighting),  then  H3  in  (141)  reduces  to  H,,  in  (133).  Alternatively,  if 
«o  ->  0  (no  frequency  modulation  in  the  waveform  s),  then  eo(q  -  1)  = 
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WW-ftK  njf.  W  UFAXW  \>  \F  rJ*  "jr  Tjr  -j/  ;rar*  ir^  vrywir*  *-  it*  v*  v*  it*  v*  ir*  vr*  yw* 


.*>i  -’V  .'v  ••  "W  ‘  ■ 


1  2  4 

1^3  ~  r  f  r  f  r  “  °o  = 


fc  •=;]•- 


°  °o> 


This  is  minimized  by  choosing  a  =  0  and  a  =  oQ,  giving  value  4 

as  usual.  Finally,  for  given  q,  H3  in  (141)  is  minimized  by 

2  2-1/2 

choosing  r  =  (1  +  (q  -  1)  e  )  ,  for  which  the  minimum  H„  = 

0  W 

2  2  1/2 

2  i-  2(1  -t-(q  -  1)  8q)  ;  however,  again,  this  increases  as 

2  fq  -  lie  as  e  increases. 

I  »o  o 


(144) 


SMOOTHING  THE  WDF 


It  was  demonstrated  in  (106)  that  the  double  convolution  of  any  two  WDFs 
is  always  nonnegative,  and  is  in  fact  equal  to  the  S1SE  of  one  waveform 
relative  to  the  other: 

tf  ii? 

Ws(t,f)  ®  Wu(t.f)  =  |Su(t,f)p  .  (145) 

This  suggests  that  one  should  choose  a  (legal)  WDF  for  weighting  u  which  is 
as  narrow  as  possible  (least  area  or  spread)  in  the  t.f  plane,  in  order  to 
minimize  the  inherent  spreading  that  (145)  implies.  The  simple  examples  in 
the  previous  section  demonstrated  that,  for  the  best  choices  of  duration  and 
linear  frequency  modulation  parameters  in  the  Gaussian  weighting,  an 
increase  of  .5  in  the  effective  area  in  the  t,f  plane  of  the  S1SE,  relative 
to  the  WDF,  resulted. 

PHILOSOPHY  AND  APPROACH 

Since  fine  detail  of  the  WDF  Ws(t,f)  will  likely  vary  in  different 
portions  of  the  t,f  plane,  this  suggests  the  following  possible  procedure 
for  analysis:  For  a  given  waveform  s(t),  compute  and  plot  the  WDF  Ws(t,f) 
according  to  (3)  or  (61).  Locate  a  t,f  region  of  interest  in  the  plane, 
where  large  (perhaps  oscillatory)  values  of  W$  occur;  denote  the  center  of 
the  region  as  tc,fc.  Estimate  the  duration,  oc ,  and  linear  frequency 


modulation  index,  ac>  of  this  particular  region  in  the  t,f  plane.  Perform 
the  S1SE  of  waveform  s(t)  according  to  ( 1 02 ) -( 1 03) ,  with  weighting 


u ( t )  =  (w  oc)  exp  - 


t2  .  ac  ,2]  2 

'  2  +  1  2  1  ‘  6c  =  ac  °c  * 

L  2cc  J 


(146) 


(for  reasons  to  be  given  below),  but  only  for  locations  t,f  in  the  plane 

near  t  =  t  ,  f  =  f  . 

c  c 


The  WDF  of  weighting  (146)  is  (with  <*c  =  2irBc) 


Wu(t,f)  =  2  exp  -  ~2  - 


t2  2  2  2 

-  ^  -  4/o2(f  -  Bct)2 

O 

L  c 


=  2  exp  -  (1  +  e2)  +  4*fte  -  4w2f2o' 

2  c  c  c 

L  °c 


(147) 


which  has  a  contour  ellipse,  at  the  1/e  relative  level,  of  area  1/2  in  the 
t.f  plane,  regardless  of  a and  «c.  This  STSE  procedure  is  equivalent 
to  smoothing  the  WDF  W^  of  waveform  s  with  the  WDF  in  (147),  for  values 
near  t  ,f  in  the  t.f  plane.  Thus  we  have  two  alternative  procedures 
for  conducting  the  smoothing  of  a  calculated  WDF  W$,  the  first  via  direct 
evaluation  of  double  convolution  (145)  for  values  of  t,f  near  regions  of 
interest,  and  the  second  via  the  STSE  in  (102)  and  (103).  Which  one  to 
adopt  will  likely  depend  on  the  number  of  points  that  must  be  closely 
investigated  in  the  t.f  plane. 


For  other  regions  of  interest  in  the  t,f  plane  of  the  original  WDF  W$ , 
different  values  of  tc,  fc>  oc>  a c  must  be  extracted  and  the  smoothing 
procedure  repeated.  Although  tedious,  this  procedure  will  minimally  spread 
the  WDF  W$  (by  area  .5)  and  it  will  guarantee  a  nonnegative  distribution. 
This  procedure  is  similar  to  that  given  in  [2];  however,  the  information 
required  to  implement  [2]  is  not  easily  available,  and  the  current  approach 
is  not  limited  to  constant-magnitude  waveforms.  A  fine-grained  analysis  of 
a  given  general  waveform  s,  for  various  t,f  values  and  yielding  nonnegative 
distribution  values,  is  not  going  to  be  achieved  without  the  expenditure  of 
considerable  effort  and  interaction  between  a  user  and  preliminary  analysis 
results . 

This  two-stage  procedure,  of  observing  the  raw  WDF  and  then  computing 
different  smoother  versions  in  different  regions,  avoids  the  arbitrary 
pre-selection  of  time  duration  and  frequency  modulation  content  of  the 
weighting  in  the  STSE,  which  would  overly  smear  the  modified  WDF  for 
improper  matches  of  parameter  values.  It  also  guarantees  nonnegative 
estimates.  In  trade,  there  is  approximately  an  increase  of  .5  in  the 
effective  area  of  the  distribution  In  the  t,f  plane  that  must  be  accepted, 
in  addition  to  a  decreased  peak  value.  For  WDFs,  W^,  with  lobes  which 
already  occupy  portions  of  the  t,f  plane  with  areas  significantly  greater 
than  .5,  this  additional  spreading  (by  area  .5)  is  not  very  damaging, 
provided  that  a and  ac  are  chosen  correctly.  Perhaps  simultaneous 
plots  of  WDF  Ws(t,f)  and  STSE  |su(t,f)|  2  would  yield  maximum 
information  about  waveform  s. 
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In  actual  practice,  where  the  integral  definition  in  (102)  is  replaced 
by  a  numerical  summation  of  samples  taken  at  increment  A,  the  quantity  (145) 
is  necessarily  approximated.  This  problem  is  addressed  in  appendix  E,  where 
it  is  shown  that  the  dominant  term  in  the  numerical  approach  is 
approximately  the  desired  quantity  (145).  Furthermore,  since  the  definition 
in  (E-l)  involves  a  magnitude-square,  the  approximation  is  guaranteed  to  be 
nonnegative.  This  need  not  be  the  case  if  the  double  convolution  of  WDFs 
and  Wu  in  (145)  is  approximated  by  sampling  directly  in  the  t,f  plane 
and  performing  a  double  summation.  However,  for  small  enough  increments  in 
both  t  and  f,  this  nonnegative  aspect  should  be  small  and  probably 
negligible;  this  latter  approach  was  used  in  [1],  although  the  smoothing 
function  was  not  a  legal  WOF. 

ALTERNATIVE  AVERAGING  PROCEDURES 


T  *  T 

Instead  of  using  R(t,T)  =  s(t  +  ~)  s  (t  -  ~)  in  (2)  as  the  instantaneous 
correlation  at  time  t  and  separation  t,  one  could  use  a  local  average,  in 
hopes  of  improving  the  correlation  and  distribution  functions.  That  is, 
consider  correlation  definition 

R(t,T)  =  V](t)  ®  R(t,r)  =  J  dt'  V](t  -  t')  s ( t '  +  |)  s* ( t '  -  |)  ,  (148) 

where  v1  is  a  fairly  sharp,  even,  real  function  centered  at  the  origin. 

The  corresponding  "locally  averaged"  WDF  is 
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W(t,f)  =  J  dx  exp( -i 2*fr)  R(t,f) 


=  fdt*  v-j (t  -  V)  w(f,f)  =  Vl(t)iw(t,f)  . 


This  is  a  convolution,  in  time  only,  of  the  WDF  of  s  with  weighting  v-j 


Reference  to  the  discussion  following  (83)  reveals  that  this  form  of 


averaging  is  inadequate,  since  it  does  not  average  additionally  on 


frequency.  Also,  (149)  need  not  remain  positive,  as  would  be  desired  of  a 


smoothed  WDf . 


Furthermore,  the  Fourier  transform  in  (149)  (as  well  as  (3))  is  over  a_H 


T,  thereby  involving  argument  values  of  waveform  s  in  (148)  which  are  very 


distant  from  the  time  point,  t,  of  interest.  If  W(t,f)  or  W(t,f)  is  to  be 


considered  as  the  "spectrum  at  time  t,"  it  is  hard  to  justify  why 


arbitrarily  distant  time  points  from  location  t  should  enter  into  their 


evaluations.  Therefore,  in  addition  to  the  local  average  in  (148)  for 


stability  purposes,  there  should  be  a  weighting  inf  in  (149)  to  better 


confine  the  Fourier  transform  to  local  values  of  waveform  s  about  time 


instant  t  of  interest. 


To  this  aim,  consider  the  more  general  form  of  average  given  by 


R(t,T)  =  v2(t,T)®  R(t,T)  =  Jdt’  v2(t  -  t'.r)  R(t\T)  = 


dt '  v,(t  -  t 


.*)  f 


df  '  exp(  i  2-wX.f  1 )  W(t'  ,f ')  , 


(150) 


W 


where  weighting  V2  depends  additionally  on  T.  Define  its  transform 

V2(t,f)  -  J  dt  exp(-i2*ft)  v2(t,t)  .  (151) 

A 

Ihen  the  modified  WDf  corresponding  to  R  in  (150)  is 

W(t,f)  --  j*|dt'  df '  W(t'.f')  V?(t  -  t\  f  -  f‘)  - 

tf 

-  W( t , f )  S>  V2( t.f)  .  (152) 

which  is  a  double  convolution  of  W  with  V2,  on  both  t  and  f.  However, 

A 

since  V2  need  not  be  a  WDF,  W  in  (152)  can  become  negative  for  some  t,f 
values.  This  form  of  smoothing  was  considered  previously  in  [9;  (1.5)]  and 
[10;  (2.1)]. 

An  additional  justification  of  two-dimensional  smoothing,  from  the 
frequency  domain  alternative  viewpoint,  is  given  in  appendix  F.  Also,  a 
general i zat ion  of  the  Gaussian  WDF  (147),  with  arbitrary  area  and  linear 

A 

frequency  modulation  content  and  which  guarantees  a  positive  distribution  W, 
is  given  in  ( F  -7) -( F -1 9 ) ;  this  result  generalizes  that  in  [11]  for  no 
frequency  modulation. 

If  we  specialize  weighting  v?  in  (150)  to  the  form 

v2(t,X)  =  u ( t  +  f)  u*(t  -  |)  ,  (153) 


then  (151)  yields 
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and  the  general  result  in  (152)  specializes  to  (106),  which  is  guaranteed 
positive.  Thus  the  special  case  of  weighting  Vp  in  (153)  leads  to  the 
STSE  of  s,  relative  to  u. 

EFFICIENT  CALCULATION  OF  SHORT-TERM  SPtCIRAL  ESTIMATE 

If  we  employ  the  weighting  u  in  (146)  with  linear  frequency  modulation 
parameter  ac,  the  spectrum  in  (102)  becomes 

Su(t,f)  *  J*  dt]  exp(-i2»ft])  s(t1)  u*(t  -  t] )  = 


=  exp(  -i2*f t) 


exp(-i2*ftp)  s(t  *  tp)  u  (-  tp)  = 


(it  a\) 


-1/4 


exp(-i2*ft)  j"  dt?  exp( -i2irf tp)  s(t 


♦  V 


expf- 


-  -\ 


(155) 


The  exp  ^-tp/(2o^)^  term  gates  out  the  portion  of  s(t  +  tp)  near 

2 

the  origin  in  tp,  while  the  exp(-i  «ct  /2)  term  cancels  linear 
frequency  modulation  in  waveform  s. 


An  approximation  to  (155)  is  obtained  by  sampling  at  increment  A  and 
using  the  Trapezoidal  rule: 
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f 

! 


-1/4 


Sft.f)  ^  (*  a)  exp(-i2*ft)  A  exp( -i2*f Ak)’ 

U  C 


*  s(t  4-  kA)  exp 


-  2  ‘*4  *  '  *c 

a 

C 


(156) 


which  has  period  1/A  in  f.  In  particular,  the  approximation  to  the  STSE, 
at  selected  points,  is 


I  2  2 

-^7  £wM2‘nk/N) 


r 


*  s(mA  +  kA)  exp 


I  1,2.2m  .  x 
"  2  k  A  fa  *  i  «Ci 


(157) 


which  is  an  N-point  discrete  Fourier  transform;  m,  n,  N  are  integers. 


The  procedure  for  analysis  is  as  follows:  for  a  region  of  interest 


centered  at  t  ,f  in  the  t.f  plane,  choose  time  values  mA  near  t  . 

c  c  c 


Then  for  each  m,  sweep  out  n  such  that  frequency  n/(NA)  is  near  f  ;  an  FFT 


will  give  all  f  values  in  (0,1/A).  Plots  of  (157)  give  a  fine-tuned  STSE 


near  t  ,f  for  the  particular  choices  of  a  ,a  .  Additional  estimates 
c  c  c  c 


with  different  parameters  will  be  required  in  other  regions;  there  is  no 
globally  optimum  smoothing  that  will  yield  high-quality  positive  spectral 
estimates  for  all  t,f  values. 


3 


The  numerical  evaluation  of  the  exponential  quantities 


Q2(k)  =  exp 


1  ,2.2/  1 

-  2  k  A  [~2  *  1QCJ 
^  \°c 


=  Q2(-k) 


(168) 


in  (157)  can  be  effected  very  efficiently  by  the  methods  given  in  [12]. 

Ihey  are  given  by  recurrences  (which  need  to  be  evaluated  only  once  for  each 


VV 


Q1  (k)  -  Q^k  -  1)  exp(2c2) 


Q2(k)  -  02(k  -  1)  Q, (k) 


for  k  >  1  , 


(159) 


with  starting  values 


0] (0)  -  exp( -  c 2 )  .  02(0)  =  1 


c  ,  -  1  ^ 
C2  2  0 


(160) 


Only  two  complex  multiplications  per  stage  are  required  in  (159). 
Furthermore,  since 


exp(-  c 2 )  =  exp 


E  (C  +  iS)  , 


(161) 


and 


exp(2c?)  - 


2  2 

C  -  S  -  12SC 


(162) 


only  one  exp,  cos,  and  sin  must  be  evaluated  to  accomplish  (159)  for  all  k. 


W0(-  WITH  MINIMUM  SPRtAD 


The  virtues  of  smoothing  WOF  W  of  waveform  s  with  the  WDF  W  , 

s  u 

(147),  of  weighting  u,  (146),  were  discussed  earlier  in  this  section.  At 
that  time,  the  selection  of  form  (146)  for  the  weighting  was  seemingly 
arbitrary.  However,  it  is  shown  in  appendix  6  that  the  weighting,  u,  which 
has  a  minimally  spread  WOF,  is  precisely  that  given  in  (146).  The  measure 
of  spread  is 


I 


II 


dt  df 


wu(t,f) 


( f  - 


Bct)‘ 


(163) 


where  (3c  is  a  specified  (observed)  slope  of  interest  in  the  t,f  plane,  and 

a  =  2*0  .  This  measure  of  spread  concentrates  the  WDF  about  the 

c  c 

specified  slope;  see  (147).  The  actual  minimum  value  of  spread  (163)  is 
given  in  ( G -24 )  as 


minimum  I  =  — jj  ,  (164A) 

8w  a 

c 


when  weighting  u  is  constrained  to  have  mean  square  duration 


[dt  t2  |u(t)|  2  =  -f  . 


( 1 64B) 


in  addition  to  unit  energy.  Without  these  two  constraints,  the  minimization 
of  spread  (163)  is  ill-posed;  see  appendix  G. 


PERFORHANCt  IN  NOISE 


In  this  section,  we  investigate  the  bias  and  stability  of  a  WDF  estimate 
obtained  from  a  noisy  waveform.  In  particular,  the  given  waveform  x  is 

x(t)  =  s(t)  +  n(t)  ,  (16b) 

where  s  is  a  deterministic  signal  of  interest,  and  n  is  an  additive 
zero-mean  stationary  noise.  In  fact,  we  have 

n(tT  =  0  .  n(t1)n*(t2)  =  Cn(t]  -  t?)  , 

Gn(f)  =  (  dT  exp(-i2irfD  Cn(T)  .  (16b) 

where  Cn  and  Gn  are  the  noise  covariance  and  power  density  spectrum, 
respectively. 

WAVEFORM  WEIGHTING 

If  the  WOF  of  given  waveform  x  In  (165)  were  directly  evaluated  via 
definition  (3),  the  result  would  be  infinite,  since  the  NxN  (noise 
cross-noise)  terms  do  not  decay  for  large  arguments.  Also,  since  the  signal 
s  will  be  assumed  to  be  transient  and  decay  to  zero  for  large  arguments, 
some  gating  or  weighting  of  given  waveform  x  is  appropriate,  in  order  to 


concentrate  on  the  time  regions  where  signal  s  is  largest.  Accordingly,  we 
consider  the  weighted  waveform 


y(t)  =  v(t)  x(t)  =  v(t)  [s( t)  +  n(t) ]  ,  (167) 

where  v(t)  is  a  deterministic  function  under  our  control. 


The  WOF  that  will  be  calculated  is  therefore 


wyy(t.f)  =  J  dt  exp( -i2«fx)  y(t  +  f)  y*(t  -  \)  = 
=  a+  b+  c+  d  , 

where 

a  =  JdT  exp( -i2irfT)  Rvv(t,x)  s(t  +  |)  s*(t  -  f)  , 

b  =  jdT  exp( -i2»fx)  Ryv(t,x)  n(t  +  |)  n*(t  -  |)  , 

c  =  JdT  exp( -i2*fT)  Ryv(t,x)  s(t  +  |)  n*(t  -  |)  , 

d  =  JdT  exp(-i2irfx)  Rvv(t,r)  n(t  +  f)  s*(t  -  |)  , 

and 

Ryv(t.x)  =  v(t  +  f)  v*(t  -  ])  . 


(168) 


(169) 


(170) 


The  first  two  quantities  in  (168)  are,  respectively,  SxS  and  NxN  terms, 

while  the  last  two  are  SXN  terms;  here,  S  denotes  signal,  while  N  denotes 

noise.  The  SxS  term,  a,  in  (169)  is  real  and  non-random,  while  NxN  term,  b, 

is  real  and  random.  On  the  other  hand,  the  SxN  terms,  c  and  d,  are  complex 
* 

random,  with  d  =  c. 


-tvjkjwv  w  wv.vjwsfww*  r*  .->  v  »ji  v  v’j.^v’j'1'  ,v  'r*  l^v 
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T^'  -  T-  ’  1 


\^YTX*^jrijrrji~j  .V-VJV'J'vT.^^rvw^ 


MEAN  VALUES 


An  alternative  expression  for  the  SxS  term  in  (169)  is 


a  -  jdt  exp( -i2.fr)  Rv>(t,T)  « ss(t ,r)  - 


-  Wvv(t,f)«Hss(t.f)  . 


in  terms  of  the  WDFs  of  weighting  v  and  signal  s. 


The  mean  of  the  NxN  term  is 


(171) 


b  -  Jdr  exp( -i2irft)  Ryv,(t,r)  CR(r)  = 


Wvy(t.f)  ®Gn(f)  , 


(172) 


where  we  made  use  of  (166).  And  since  noise  n  has  zero  mean,  there  follows, 
for  the  SxN  terms,  ~  -  <f  =  0.  Collecting  these  results  together,  the  mean 
of  WDF  estimate  (168)  is 


Wyy(t.f)  =  Ww(t*f)e[Wss(t’f)  +  Gn(f)] 


(173) 


No  additional  statistical  properties  on  the  noise  n,  such  as  a  Gaussian 
process,  are  needed  for  result  (173);  this  holds  for  an  arbitrary  stationary 
noise  process.  The  difference  in  mean  outputs,  for  signal  present  versus 
signal  absent,  is  just  a,  as  given  by  (171). 


VAR L ANC t  OF  WOF  ESTIMATE 


In  order  to  determine  the  variance  of  estimate  (168) ,  we  need  to  assume 
that  noise  n  is  Gaussian.  Furthermore,  in  addition  to  properties  (T66),  we 
will  presume  that 

n( t1 )  n(t2)  =  0  ,  (174) 

as  is  true  when  n  is  an  analytic  process  or  a  complex  envelope  [13,  ch .  2  ] . 
Then  (168)  yields 

Wyy(t,f)  =  |wyy(t,f)l  2  =  +  b2  +  lcl2  *  )dJ2  f 

+■  2  ab  +  cd*  +  c*d  ,  (175) 

the  other  terms  being  zero  due  to  n  being  zero-mean  Gaussian  noise. 


The  second  term  on  the  right-hand  side  of  (175)  can  be  developed  from 
(169)  as 

b2  =  Ibl 2  =  J[dT1  dt2  exp  (-12*^,  -I2j)  Rvv(t,T1)* 


*Rvv^*T2MCn(ri) 


VT1 


C*(TJ  +  C2 


(176) 


where  we  used  (166)  and  (174).  Referring  to  (172),  we  have 
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i*i.i** 


b2  -  b  2  +  J/dr1  dt2  exp  (-12^^  -Tp)  Rvv(t,T1)’ 


**„l*'\)  c 


At  this  point,  it  is  convenient  to  define 


Gn2)(f)  =  JdT  exP( -i4irfc)  C2(I) 


Gf,2)(f/2)  =  JdX  exp(-12*fT)  C2(X)  =  Gn(f)«>Gn(f) 


C2(X)  =  J  df  exp(i2trfT)  6^2>(f/2)  = 


2  J  dw  exp(i4irwT)  Gj,  '(v) 


When  this  result  is  substituted  In  (177),  there  follows 


-  ^ 

b2  =  F  +  2  Jdv  W2y(t,f  -  v)  G^2) ( v)  = 


2  f 

=  r  +  2  W2v(t,f)®  G<2)(f)  . 


v k\i  .  ■ v  v  *  .*'*  *  v** 


k 
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The  third  and  fourth  terms  in  (175)  are 


|di 2  =•  ]c|2-  {{<*]  ^2  exp(-i2trf(T1  -T2))  Ryv(t ,1^)’ 


where 


*  Ryvtt.fj)  S 

=-  {  dv  Gn(v)  J  B^t.f  -  , 

B(t.f )  2  J  dZ  exp(-i2trfr)  Rvv(t,r)  s  (t  +  = 

=  J  dv  exp(  i 2rirt)  Wyv(t,f  -  ^  S(v)  . 


(182) 


(183) 


(As  special  cases,  if  weighting  v(t)  =•  1  for  all  t,  then  Ryv(t,T)  =  1 


and 


|c| 2  =  4  {dv  Gn(v)  |s(2f  -  v)|  2  ; 


while,  instead,  if  Gn(f)  -  for  all  f,  then  Cn(T)  =  6(r)  and 

\2 


|c|2  =  2  Nd  J  a  R2vv(t,T)  |s(t  *  || 


If  both  conditions  above  hold,  then 


(184) 


(185) 


yj 

M 


|c|  2  =  4  N,  E  .  ) 


(186) 
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Returning  to  the  general  case  again,  the  fifth  term  in  (175)  is  given  by 


combining  (171)  and  (172),  while  the  sixth  term  is 


cd*  =  ||  dt1  dr2  exp(-i2irf(T1  -T^))  R^t,^) 


*  Rvv(t.T2)  s 


^  s*^  -  ^  t\*(t  -  3j^n*^t  f  =  0  ,  (187) 


by  use  of  (174) 


Combining  the  above  results,  we  have,  for  the  variance  of  the  WDP 
estimate, 

Var{wyy(t,f)|  =  a2  *  b2  +  2  |c| 2  ♦  2ab  -  (a  +  b)2  - 

__  _2  _ 

=  b2  -  b  +■  2  1c  j 2  = 

-  2  wjv(t,f)  6  G<2)(f)  +  (NxN) 

2 

+  2  |  dw  Gn(v)  |B(t,f-f)|  .  (SxN)  (188) 

This  result  holds  for  arbitrary  signal  s,  weighting  v,  and  noise  spectrum 
Gn.  The  quantities  G^  and  B  are  defined  in  (178)  and  (183), 
respectively. 
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If  we  do  not  weight  waveform  x(t),  that  is,  choose  v(t)  =  1  for  all  t  in 

(167),  then  R  (t,T)  =  1,  Wvv(t,f)  =  6(f),  and  the  NxN  term  in  (188) 

becomes  infinite;  that  is,  the  WDF  estimate  (168)  has  infinite  variance  if 

we  do  not  weight  in  time,  regardless  of  what  the  actual  noise  spectrum, 

G  ,  is. 
n 


On  the  other  hand,  if  the  noise  n  is  white,  then  Gn(f)  -  N^  for  all  f, 

C  (Z)  =  N.  6(T) ,  and  in  (178)  is  infinite,  which  makes  the  NxN  term 

n  d  n 

in  (188)  infinite.  Thus,  if  we  do  not  filter  out  the  noise  which  is  out  of 
the  band  of  the  signal,  the  WDF  estimate  has  infinite  variance,  regardless 
of  what  time  weighting  v  is  employed. 


WDF  PROCtSSOR 

In  view  of  the  above  observations,  we  now  consider  the  general  WDF 
processor  depicted  in  figure  7.  The  only  new  element  here  is  the 


Figure  7.  WDF  Processor 


t ime -invariant  linear  filter  with  transfer  function  H.  The  input  noise  n( 
is  presumed  to  be  white  over  the  band  of  the  input  signal  sQ;  mathematically, 
this  is  handled  by  letting 
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where  is  the  double-sided  noise  spectral  density  level  in  watts/Hz. 

The  linear  filter  H ( f )  approximately  matches  the  bandwidth  of  the  input 
signal  and  passes  SQ(f)  essentially  unaltered,  while  filtering  out 
undesired  noise  spectral  components.  The  actual  filter  output  signal  s  is 
given  by 

s(t)  =  h(t)®  s0(t)  =  Jdf  exp(12*ft)  H(f)  SQ(f)  .  (190) 

The  weighting  v(t)  approximately  matches  the  duration  of  the  signal  and 
passes  s(t)  essentially  unaltered,  while  gating  out  undesired  noise  temporal 
components.  Representative  plots  of  the  various  quantities  in  figure  7  are 
given  in  figure  8. 

A  numerical  example  of  the  WDF  processor  in  figure  7  is  carried  out  in 
complete  detail  in  appendix  H,  including  the  mean  and  variance  results  given 
earlier  in  this  section.  In  particular,  the  input  signal  sQ  is  a  linear 
frequency  modulation  waveform  with  Gaussian  amplitude  modulation. 


SUMMARY 


When  a  segment  of  a  stationary  random  process  is  available,  the  method 
of  Blackman  and  Tukey  [14]  tells  us  that,  to  estimate  the  correlation 
function  at  delay T,  we  should  average  the  product  of  waveform  values 
separated  in  time  by  X  seconds,  and  that  we  should  carry  out  this  averaging 
over  the  total  available  data  record,  in  order  to  reduce  the  effect  of 
random  fluctuations.  For  a  nonstationary  process,  the  averaging  interval  is 
further  limited  to  that  in  which  a  significant  change  is  statistics  does  not 
occur. 

After  obtaining  the  estimated  correlation,  the  Blackman-Tukey  method 
further  directs  us  to  weight  the  correlation  values  in  the  neighborhood  of 
T=  0  more  heavily  than  those  for  larger  T,  and  to  Fourier  transform  the 
weighted  correlation  estimate.  The  weighting  should  taper  off  to  zero  for 
larger!,  so  as  to  suppress  these  more  noisy  estimates,  and  the  taper  should 
be  gradual  so  as  not  to  create  significant  positive  and  negative  sidelobes 
in  the  frequency  domain. 

These  two  operations,  averaging  in  time  and  weighting  in  delay,  are  both 
totally  absent  in  the  WOF,  as  may  be  seen  from  (2)  and  (3).  In  fact,  (2) 
and  (3)  might  be  viewed  as  the  ultimate  in  greediness  of  a  spectral 
estimate,  since  they  include  no  averaging  and  no  weighting.  Viewed  in  this 
light,  It  Is  not  surprising  that  the  WOF  has  some  very  debilitating  behavior 
in  terms  of  negative  distribution  values  and  large  Interference  terms. 


a*A  «*i 


»,«  i 
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The  inclusion  of  averaging  and  weighting  in  the  spectral  estimate,  as 
typified  by  ( 1 50) -(152)  and  (F-l),  results  in  a  modified  distribution 
function  which  is  a  double  convolution  with  a  smoothing  function  in  the  t,f 
plane.  Furthermore,  the  averaging  and  weighting  in  (F-l)  takes  place  both  in 
the  f  and  v  domain  (line  4)  just  as  well  as  in  the  t  and  X  domain  (line  3). 
Alternatively,  line  2  indicates  that  the  complex  ambiguity  function  may  be 
weighted  in  two  dimensions  and  doubly  Fourier  transformed.  However,  the 
resultant  modified  WDF  need  not  be  positive. 


The  identity  of  this  double  convolution  with  a  positive  STSE,  when  the 
smoothing  function  is  a  legal  WDF,  allows  for  an  alternative  approach  that 
is  very  attractive  computationally  and  is  easy  to  interpret.  The  preliminary 
calculation  of  the  WDF  serves  to  point  out  regions  of  interest  in  the  t,f 
plane  and  to  quantify  the  time  and  frequency  extents,  as  well  as  the  amount 

of  linear  frequency  modulation,  to  utilize  in  weighting  u  in  the  STSE.  This 

.  .  2 
procedure  is  illustrated  in  appendix  I  for  the  waveform  s(t)  =  t  exp(-t  /2) 

and  shown  to  yield  a  physically  meaningful  smoothed  distribution  function, 

whereas  the  WDF  is  very  difficult  to  justify  and  interpret  on  any  physical 

grounds . 


It  was  pointed  out  earlier  that  double  convolution  of  a  given  WDF  with  a 
Gaussian  WOF  increases  the  spread  of  the  smoothed  function  by  area  .5  in  the 
t,f  plane,  since  the  effective  area  of  a  Gaussian  WDF  is  .5.  Strictly 
speaking,  this  is  only  true  when  the  Gaussian  WDF  contour  ellipse  has  the 
same  tilt  and  the  same  ratio  of  major-to-minor  axes  as  the  given  WDF 
(assumed  Gaussian  in  the  region  of  interest  in  the  t,f  plane).  More 


generally,  if  there  is  a  mismatch  in  tilt  or  ratio  of  major-to-minor  axes, 
the  effective  area  is  increased  by  more  than  .5,  thereby  leading  to 
additional  spreading  in  the  t,f  plane.  The  detailed  derivations  are 
presented  in  appendix  J. 

The  performance  of  an  estimator  of  the  WDF  of  a  signal  in  the  presence 
of  noise  depends  on  the  amount  of  filtering  and  weighting  employed  to 
suppress  noise  components  in  frequency  and  time.  Exact  relations  for  the 
mean  output,  the  bias,  and  the  variance  of  the  WDF  estimate  are  given. 


TR  8225 


Ja <  I. 


APPENDIX  A.  SLICES  IN  TIME  OF  THE  WDF 


The  voltage  density  spectrum  of  waveform  $(t)  was  given  in  (9).  If  s(t) 
is  sampled  at  increment  A^.  an  approximation  is  afforded  according  to 

S{f)  =  j"  dt  exp(  i2wft)  s(t)  = 

S  At  exp(  ~i2irfkV  s(kAt)  -  $(f)  for  all  f.  (A-T) 


The  summation  on  k  runs  over  the  range  of  nonzero  summand.  Since 


H) 


-  S( f )  . 


(A-2) 


then  5(f)  has  period  1/At  in  f.  We  limit  the  evaluation  of  S(f)  to  the 


values 


2 


exp( -i2«nk/Nfr)  s(kA^)  for  0  <  n  <  N^  -  1  ,  ( A-3 ) 


where  n  and  Nf  are  integers,  and  thereby  cover  a  full  period  of  S(f).  A 
representative  plot  of  )s(f)l  and  its  sampled  values  appears  in  figure  A-l 


For  the  low-pass  case  of  5(f)  depicted  in  figure  A-l,  it  is  necessary  to 


choose 


At  <  (2fH)  1 


(A-4) 


1  *vv".y  T'  *'*[*%** 


.s  V  -‘4-4  */  j.  '%■  .. .  ?.'  1' <  ■ 


TR  8225 


(V.) 


0  i  Kf 

Figure  A-l .  Low  Pass  Spectrum  ^(f) 


in  order  to  avoid  aliasing.  We  will  also  need  frequency  spacing 


(NfAt)  1  <  (2T)'1  ,  that  is,  Nf  >  2T/At 


(A-5) 


in  order  to  track  the  wiggles  of  $(f)  in  frequency,  where  1  is  the  effective 
time  duration  of  s(t).  In  fact,  we  may  need  frequency  spacing  (NfAt)  1  to 
be  very  small  if  we  are  to  do  further  accurate  integrations  on  S(f).  Thus 


we  need 


Nf  >  4  fH  T, 


( A-6) 


and  perhaps  much  larger  for  further  manipulations. 


The  WOF  can  be  written  in  terms  of  the  spectrum  S(f)  according  to  (10): 


[t.f)  =  j  dv  exp(i2irvt)  S(f  4-  |)  S*(f  -  f)  = 


-I 


du  exp(i4irut)  S(f  +  u)  S  (f  -  u)  . 


(A-7) 


.V 


2 

§ 


IftV 

ft 

is 


I 
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If  we  sample  S(f)  at  increment  Af  in  frequency,  the  Trapezoidal 
approximation  to  the  WDF  is 

W(t.f)  =  2Af  2.  exp(  i4«tJ?Af )  S(f  +^Af)  S*(  f  -X  Af) 


( A-8) 


for  all  t.f .  Since 


-  W(t.f)  , 

A  .] 

then  W(t,f)  has  period  (2Af)  in  t.  Accordingly,  we  evaluate  only 
W^flf  .  nAf^  -  2Af  ^exp(i2»mi/Nt)  s((n  )Af)  S*((n  -J?  )Af) 


( A -9 ) 


( A -1 0) 


for  0  <  m  <  Nt  -  1,  where  m,  n,  Nt  are  integers.  In  this  manner,  we  get 

A 

a  slice  of  W(t,f)  in  time  t  (m)  for  fixed  frequency  f  (n).  The  operation  in 
( A  - 1 0 )  can  be  efficiently  realized  as  an  N^-point  f FT  of  collapsed  samples 
when  is  a  power  of  2. 


Now  the  only  information  on  S(f)  that  we  have  available  are  the  samples 
of  S(f)  given  in  (A-3).  If  we  choose,  without  loss  of  generality. 


af  ~  (Nfat) 


(A-ll) 


then  (A-10)  becomes  (exactly) 


rt  •  v;)'  »A^-exp(12’M/" t>  s(vf)  s*(S^ 


(A-12) 


for  0  <  m  <  -  1 .  We  then  adopt  as  our  approximation  to  W,  which  itself 

is  an  approximation  to  W,  the  quantity 


k5»XX> Kyji/JSCCI  vXXX  Vf  -  IV.’*V.'iV  -  bUyXCXX'WCA. 
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“(i1  ^  '  V^) r  "A  m  exp(12,ni,/Nt)  S(M?)  5*(”,^)  (A  13) 


for  0  <  m  <  -  1 ,  where 


X  c  l.n  denotes  |n  +jt|  <  Nf/2  , 


(A -14) 


and  we  presume  that  S ( f )  has  been  calculated  for  |f|  <  (2a^)  ^ ;  see 


figure  A-l 


The  quantity  may  have  to  be  large,  in  order  to  sample  S(f) 


finely  enough  for  an  accurate  WDF .  Then  since  the  spacing  in  t,  applied  to 


Nf  Af 

W(t,f),  is  - .  it  may  require  a  large  value  of  in  order  to  keep  track 


of  the  variations  versus  t.  Also,  n  may  not  have  to  run  through  consecutive 


integer  values,  but  may  take  on  decimated  values,  so  that  n/fN^) 


tracks  the  f  behavior  adequately. 
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APPENDIX  B.  OSCILLATING  WDF  FOR  SEPARATED  PULSES 


Consider  the  waveform  s(t)  in  figure  B-l  consisting  of  two  separated 


energy  bursts  of  general  shape.  Let  t^  and  t?  represent  the  "center"  of 


each  pulse,  and  let  T1  and  be  some  measure  of  their  durations.  Define 


tc  =  2^1  f  ^2^  ' 


(B-l) 


Let  us  investigate  the  WDF  of  s(t)  for  t  near  tc>  that  is,  near  the  center 


of  the  two  pulses.  In  particular,  let  time 


t  =  tc  +  A  , 


(B-2) 


where  a  is  small.  Then  from  (3), 


W(tc  *  A , f )  =  J  dl  exp( -12«ft)  s(tc  +  A  *  *)  s*(tc  A  -  1> 


r  exp  (-i2*f (t-j  -  t2))  J  dx  exp(  -i2*fx)  s(t2  +  A  +  |)  s*(t1  +  a  -  |) ,  (B-3) 


Figure  B-l.  Waveform  s(t) 


4i  V.\*  ,«v>v\,.v 
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where  we  let  T=  +  x. 

It  can  be  seen  from  the  integral  expression  in  (B-3)  and  figure  B-l , 
that  for  small  A,  only  small  values  of  x  will  contribute  to  the  value  of  the 
WDF .  In  fact. 

|x|  <  min(I1,72)  (B-4) 

is  the  dominant  range  of  contribution  to  the  Integral.  Thus  the  variation 
with  f  of  the  Integral  component  of  (B-3)  Is  slow.  By  contrast,  the  leading 
exponential  term  In  (B-3)  varies  much  faster  with  f,  since 

t2  -  t!  >  max(T1.T2)  >  Rln(T1(T2)  .  (B-5) 

Since  these  faster  varying  oscillations  of  the  exponential  term  cannot  be 
cancelled  by  the  slower  integral  contribution  of  (B-3),  the  WDF  will 
oscillate  in  f,  for  times  t  near  tc  *  (t^  *■  t2)/2.  Thus  separated 
time  pulses  will  lead  to  oscillations  (in  f)  of  the  WDF,  near  times  midway 
between  the  pulses,  regardless  of  their  detailed  shapes.  An  analogous 
argument  can  be  presented  for  spectral  components,  based  upon  form  (10)  of 
the  WDF. 

Notice  that  as  t^  approaches  t1 ,  and  the  two  pulses  become  one,  the 
oscillating  exponential  term  in  (B-3)  disappears,  allowing  for  the 
possibility  of  a  slowly  varying  (hopefully  positive)  lobe  in  the  WDF. 


8-2 
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APPtNOIX  C.  AMBIGUilY  FUNCTION  OF  (79) 


The  complex  ambiguity  function  of  waveform  s ( t )  In  (79)  Is  obtained  by 
substitution  into  (22)  and  use  of  (68): 

£(»>.T)  =  Vi  °k|  expl"  4  ^2  “  4  °K*  VKI  * 

M  L  °KI 


2  _  o2 

4-  i  j  “2"  2  1WVkt  +  l2trfkj TkS 

aX  +  °k 


-%vk»  *12’Vtk-v]  • 


(C-1) 


where 


V  • 


-2  1,  2  .  2, 

°kt  “  2(ok  *  } 


Wl 


1 


~  2 
°Kl 


2(fk  +  V  • 


'  *6  *  «i) 


Tk*  =  T  *  i  ~  H  •  vkl  2,r(v  +  i  *  V 


(C-2) 


The  diagonal  terms  of  (C-1)  are 


Ek  H-4^-i»k4'2"2tt2*fkT-12'‘k 
k  L  °k 


']• 


(C-3) 


which  are  complex  and  oscillate  with  T  and  v  due  to  the  imaginary  terms. 

The  contour  at  the  1/e  relative  level,  of  the  magnitude  of  the  k-th  term,  is 
an  ellipse  with  axes  twice  as  large  as  those  depicted  in  figure  4.  In 
addition,  the  peak  amplitude  is  decreased  by  a  factor  of  2  below  that  for 

C-1 
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the  WDF  in  (80).  Thus  the  ambiguity  function  is  a  more  smeared  function  of 
time-frequency  than  the  WDF. 


The  ambiguity  function  has  peaks  of  value 


P1  V*  Sw  «xP[i.(fk  *  i)(tk  -  t,i] 


(C~4) 


centered  at 


(T.wJ-Ufc-  Vfk-J> 


for  all  k,J?.  The  phases  of  ( C -4 )  are  virtually  random  relative  to  each 
other.  A  slice  in  v,  for  fixedT,  varies  (in  addition  to  the  Gaussian 


(0-5) 


envelope)  as 


^sbi! T**  ’ tw)  ■ 


(C-6) 


which  could  be  either  a  slow  or  fast  variation,  depending  on  the  particular 
parameter  values.  All  of  these  features  make  physical  interpretation  of  the 
ambiguity  function  very  difficult. 


\Vf\v,y,v. 
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APPENDIX  D.  ROTATION  OF  AXES 


Consider  the  general  second-order  curve  described  by 

2  2 

Ax  +  Bxy  +  Cy  +  Ox  +  Ey  *  F  =  0  . 

If  we  rotate  the  x,y  coordinate  axes  according  to  figure  D-l ,  we  have 

x  =  x'  cos(B)  -  y1  sin(B)  , 

y  =  x'  |1n(|)  f  y'  cos(B)  . 

Substitution  in  (0-1)  yields 

A'x'2  f  B'x'y'  +■  C * y ' 2  +  0 1 x 1  +  E'y’  +  F  =  0  , 


(D-l) 


(D  2) 


(D-3) 


where 


* 

\ 

\ 

1 

y 

\ 

A 

\ 

1 

^  "  V  x 

-  -  "o 

\ 

Figure  0-1.  Rotated  Coordinate  Axes 
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A'  -  ^ C A  +  C  +  (A  -  C)  cos(2B)  f  B  sin(28)] 

B'  =  B  cos(2B)  +  (C  -  A)  sin(2B) 

C'  =  j[A  +  C  -  (A  -  C)  cos(2B)  -  B  sin(2B)] 

O'  =  0  cos(B)  +  E  sin(B) 

E'  =  -  0  sin(B)  +  E  cos(B)  .  (04) 

If  we  want  to  eliminate  the  cross  product  term  In  (D-3),  we  must  make 
B*  =  0,  that  is,  take 

tan(2B)  =  •  (0-5) 

We  will  also  choose  2B  In  the  principal  value  range: 

2B  =  arc  tan  j  ^  ;  (D-6) 

that  is 

-  |  <  2B  <  |  ,  “  J  B  <  f  *  (0-7) 

All  other  solutions  for  2B  differ  by  n»;  that  Is,  B  differs  by  mr/2.  These 
are  the  major  and  minor  axes  of  the  curve  described  by  ( 0-1 ) . 

If  we  now  define 

R  =^(A  -  C)2  ♦  BZ\  P  =  sgn(A  -  C)  ,  (0  8) 

where  V"1  denotes  the  positive  square  root,  we  find 


1/  i.1  », 
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cos(2B)  =  ^  ~  Cl 


(0  9; 


see  figure  0-2.  And  since  sin(2(3)  has  the  same  polarity  as  tan(2(l)  in  the 
principal  value  range. 


sin(2B)  -  sgn^~^~^  =  |  p  • 

Also,  since 

cos2(8)  =  |(1  4-  cos ( 2B) )  =  1^  + 

then 

cos(B)  • 

sfn(B)  59"(r?-c)  ' 

It  then  follows  that  the  coefficients  in  (D-4)  simplify  to 

A'  =  1(A  4-  C  4-  RP) 

B 1  -  0 

C'  =  |(A  4  C  -  RP)  , 


(0-10) 


(0-11 ) 


(0-12) 


(0-13) 


figure  0-2.  Triangle  Interpretation  of  (0-5) 


i] 


tifvt 
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from  which  there  follows 


1  2 

AT  -  AC  -  J  B  . 


(D  14) 


Additionally,  we  have 


tan(B)  = 


cot(O)  - 


_  PR  -  (A  -  C) 


PR  4-  (A  -  C) 


(0-15) 


As  a  result  of  the  above,  the  general  equation  in  (D-3)  can  be  written  as 


:V « ail  el2 

')  4A'  4C'  h  * 


(0-16) 


where 


1  A  +  C  -  RP 


J_  -  A  +  C  +  RP 


2A'  '  4AC  -  B2  *  2C'  "  4AC  -  B2  * 


(0-17) 


The  simplest  expressions  for  O'  and  E'  appear  to  be  those  given  in  (0-4),  in 

.  .  .  2  2 
conjunction  with  (D-12).  However,  O'  and  E'  can  be  simplified. 


resulting  in  expression 


A'/x'  + 


where 


(0-18) 


~  AE2  +  CD2  -80E  _ 

G  -  p  "  F  ’ 

4AC  -  B* 


(D-19) 


Now  suppose  that  A  and  C  in  (0-1)  are  positive  and  that  4AC  >  B  . 
Then  A'  >  0  and  C'  >  0,  meaning  that  (0-18)  is  an  ellipse  if  G  >  0. 


t  »  4/  4  1  ♦  *  r 


t*  >  v\.vm 


V  i  V  /VAV  v 


»v>y  a  v  *  v  ‘Try!  • '  o  7t  r  *  ‘yys  y  , 


EXAMPLE 


Consider  the  ellipse  in  (91),  for  which 

A  =  1  +  02  ,  B  =  -  20  ,  C  ---  1;  0  >  0  .  ( D  -2 1 ) 

O  O  O  ~ 

Then  (0-5)  yields 

tan(28)  =  -  2/0o  ,  (0-22) 

from  which  there  follows 

tan  ^20  +  f )  =  -l/tan(20)  =  0q/2  , 

„  1  *  /eo\  ir  (0-23) 

8  =  ^  arc  tan  "  4  • 

As  0Q  varies  from  0  to  +«6  ,  0  varies  from  -tr/4  to  0;  thus  0  always  lies  in 
the  principal  value  range,  as  required  by  (0-7). 

However,  since  A  >  C  in  (0-21),  then  0  is  the  angle  in  the  x,y  plane  of 
the  minor  axis  of  ellipse  (91).  The  major  axis  angle  in  the  x,y  plane  is 

ir  =  r  +  \ ,  (0-24) 

which  varies  from  w/4  to  w/2.  There  follows 

tan(2^  =  tan(20  +  ir)  =  tan(20)  =  -  2/0q  (D-25) 

from  (0-22),  whereupon  the  slope  of  the  major  axis  in  the  x,y  plane  can  be 
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ta"t-  *l/9o  *  4)' 


( D -26) 


This  slope  varies  from  1  to  +«6  as  0Q  varies  from  0  to  +*>.  Conversely, 
given  a  measured  slope,  tan1^,  of  the  major  axis  of  a  WDF  contour  in  the  x,y 
plane,  the  corresponding  amount  of  linear  frequency  modulation  can  be 
determined  from  (0-25)  or  (0-26)  as 


e0  =  tan ^  -  1/tan*^"  . 


(0-27) 


The  final  determination  of  frequency  modulation  parameter  aQ  in  (84) 
requires  the  additional  knowledge  of  aQ  in  (90). 


In  practice,  where  both  oQ  and  aQ  are  unknown  apriori,  the  WOP 
will  likely  be  plotted  directly  on  the  t,f  plane.  According  to  (78)  and 
(0-24),  the  major  axis  will  then  lie  on  the  line 


f  = 


t 


tanNr 


( D-28) 


which  can  be  observed  and  measured.  But  oQ  can  be  determined  separately 
from  a  slice  in  f  (at  fixed  t)  of  the  WDF,  since  the  variation  in  f  in  (91) 
is  proportional  to 


exp 


[-  °lk'<  -  v)?] 


( D-29) 


Thus  the  distance,  between  frequency  values  that  are  down  by  1/e  from  the 
peak  on  this  frequency  slice,  is  (*o0)  1 ,  and  can  be  used  for  direct 
calculation  of  aQ.  Then  (0-28)  and  (0-27)  yield  tan^  and  eQ, 
respectively. 


0-7/0  -8 
Reverse  Blank 


to  get 


i(t.f)  =  a2  exp(  -i2i»f  Am)  R^.mA^  Ry  (t  -  n|,  mA^ 


m+n 

even 


Define,  for  use  below,  the  function 


C( t1  ,t,f )  =  ws(trf)  ®  wu(t  -  trf) 


=  j"df1  W$(t1,fI)  Wu(t  -  trf  -  f,)  , 


(E-3) 


(E-4) 


E-l 
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I 

I 

I  which  Is  the  convolution,  on  f,  of  WOfs  W  and  W  .  When  the  Integral  on 

I  s  u 

t1  Is  effected,  it  yields  the  desired  double  convolution: 


tf 

dt1  C(trt.f)  =  W$(t.f)  9  Wu(t,f)  .  (t-5) 


We  now  express  and  Ru  In  (£-3)  in  terms  of  the  inverse  transforms 
of  W$  and  Wu>  respectively,  according  to  (24),  interchange  summations 
and  integrations,  and  use  the  facts  that  (4,  chapter  2] 


^  exp(  -i2*xm)  *  ^2l4(x  -  . 

A  ii  An  ^ 


m  even 


^  exp(-12*xm)  =  ~2,(-l  f  6(x  -  ^|)  , 
riHH  X 


(E-6) 


m  odd 


to  get  approximation  (£-1)  in  the  form 


Jft.n  5.  '(“f-t.f  -&)♦ 

X  n  even 

*  -  2i)  ■  (t-V 

X  n  odd 

The  ^=0  terms  together  give,  where  the  sum  is  now  over  all  n. 


£-2 


which  Is  a  discrete  approximation  to  desired  quantity  (E-S). 


The  J?  =  1  terms  in  (E-7)  yield 

f£(-l)n  C(2|,t.f  -  ^  .  (E-9) 

II 

which  is  approximately  zero.  A  similar  result  holds  for  X  ~  *1  • 

The  Jl^  2  terms  in  (E-7)  are 

f^C^.t.f  -  \)  ,  (E-10) 

n 

which  is  a  discrete  approximation  to  convolution  (E-5),  but  shifted  by 
frequency  T/fi.  For  small  sampling  increment  &  in  approximat ion  (E-7),  the 
quantity  (E-10)  will  be  small  In  the  fundamental  region  centered  at  f  =  0, 
and  can  be  neglected.  Thus  (E-8)  Is  the  dominant  term,  giving 

J(t.f)*  ^C(n-f.t.f)  = 

n 

r  tf 

S  J  dt1  C(trt,f)  =  Ws(t,f)  6  Wu(t,f)  .  ( E-l  1 ) 


E-3/E-4 
Reverse  Blank 


^\«yitvnryvvvwi>sriA^v-\.v\^'  k~  vi  k-  wv  icro  Mqiwif’v»  »v  vj  yy  *v  y  j  w  y>  vj^.t.'?  to  wj  ru»  rgrw  mi  m  pit  **  *  w  -  u  up,  y,  irr 
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APPENDIX  F.  SOME  SMOOTHING  C0NSIDERA1 IONS 

Ibis  discussion  complements  and  extends  that  given  In  (148)-(152) 
regarding  two-dimensional  smoothing  of  the  WDF.  For  easy  reference,  we 
repeat  the  diagram  under  (24)  and  furnish  an  additional  one  for  the 
smoothing  functions  that  will  be  employed  here.  An  arrow  denotes  a  Fourier 
transform. 

WAVEFORM  FUNC1 IONS 


ft  S*(t  -  |)  =  R(t,T)  - f - W(t,f) 


f 

X(«.T)  - - A(v.f)  «  S(f  *  f)  S*(f  -  fj 


SMOOTHING  FUNCTIONS 


v2(t,r)  - r 


V2(t,f) 


<T2^v*r) 


02(v,f) 


By  using  the  basic  Fourier  transform  relations  above,  we  may  readily 
show  that  two-dimensional  smoothing  of  WDF  W  with  general  function  V?  may 
be  written  in  several  alternative  forms: 


F-l 


lTiTT 
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W(t.f)  •  V2(t.f) 


=  jjdv  dr  exp[i2ir(ut  -  fr)]  ")C{v,x)  q2(v,x)  = 


j*  dcexp(-i2irfr)  [R(t,r)  «  v2(t,x) ] 


J"dw  exp(i2«wt)  [ A( w , f )  0  Q2(v,f)] 


(F-D 


The  second  line  says  that  the  ambiguity  function  X  of  waveform  s  should  be 
weighted  by  and  the  product  then  double  Fourier-transformed  into  the 
t,f  plane.  The  third  line  indicates  smoothing  of  R  on  t,  followed  by 
transformation  of  a  weighted  function  of  x.  The  last  line  performs 
smoothing  of  A  on  f,  followed  by  transformation  of  a  weighted  function  of  «. 
These  relations  extend  those  given  in  ( 1 50) -< 1 52 ) .  The  function  V2  above 
need  not  be  a  legal  WDF.  The  volume  under  the  smoothed  distribution  (F-l) 
is  the  product  of  the  volumes  under  W  and  V2>-  if  the  latter  volume  is 
unity,  the  energy  of  waveform  s  results  again,  as  is  desired. 


INADEQUACY  OF  TIME  SMOOTHING  ALONG 


Consider  the  special  case  where  smoothing  function  V?  is  a  delta 
function  of  f;  then 

V2(t,f)  =  Vl(t)  4(f) 

v2(t,T)  =  v.|  ( t) 
q2(w.x)  =  ( v) 


02(v,f)  =  V.j  ( v)  4(f) 


r*vtf 


(F-2) 
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Equation  (F-l)  then  simplifies  to  an  averaging  of  the  WDF  solely  in  time: 

t 

W(t.f)  ®  v^t)  = 

a  fj*d»  dr  exp[i2»(wt  -  fz)]  Z(v,t)  V^w)  = 

=  J  dr-  exp(-i2*fr)  [R(t,r)  ®  v^t)]  = 

=  j*  dv  exp( i2trvt)  A( v, f )  V^v)  ,  (F-3) 

which  is  an  extension  of  ( 1 48) -( 1 49 ) . 

The  advantageous  feature  of  locally  averaging  the  instantaneous 
correlation  R  in  time,  indicated  in  line  3,  is  equivalent  to  weighting  the 
"local  spectrum" 

A(v,f)  =  S(f  +  *)  S*(f  -  *)  (F-4) 

in  line  4  by  function  V^v),  prior  to  Fourier  transforming  back  into  the  t 
domain.  This  weighting  on  v  is  sensible,  since  if  WDF  W(t,f)  or  some 
modified  version  is  to  represent  the  spectrum  at  f,  the  transform  on  »  in 
line  4  of  (F-3)  ought  not  to  involve  arbitrarily  distant  values  of  v, 
otherwise,  waveform  spectrum  S  In  (F-4)  will  then  be  utilized  at  argument 
values  very  different  from  the  frequency  f  of  interest  and  would  be 
nonrepresentative.  However,  there  is  no  weighting  onT  in  line  3  of  (F-3), 
thereby  allowing  arbitrarily  distant  argument  values  of  signal  s,  from  the 


F-3 


time  instant  t  of  interest,  to  be  considered;  this  unrealistic  feature  of 
the  WDF  is  one  of  the  reasons  for  its  undesirable  properties. 


INADEQUACY  OF  FREQUENCY  SMOOTHING  ALONE 

Now  consider  the  alternative  special  case  where  in  (F-l)  is  a  delta 
function  of  t;  then 

V2(t.f)  =  4(t)  V^f) 
v2(t,r)  =  s(t)  v^r) 
q2(«,T)  =  v^x) 

Q2(v.f)  =  V^f)  .  (F-5) 

Equation  (F-l)  then  simplifies  to  an  averaging  of  the  WOF  solely  in 
frequency: 

f 

W(t,f)  •  V^f)  = 

=  JJdw  dr  exp(i2w(wt  -  ft)]  Xi».T)  v^r)  = 

=  j"  dt  exp(-i2wft)  R(t,T)  v^r)  = 

=  |d»  exp( i2wvt)  [A(u,f)  ®  V-, (f ) ]  . 


(F-6) 


The  advantageous  feature  of  averaging  the  "local  spectrum"  A  In 
frequency,  indicated  in  line  4,  is  equivalent  to  weighting  the  instantaneous 
correlation  R(t,r)  in  line  3  by  function  v^T),  prior  to  Fourier 
transforming  back  into  the  f  domain.  This  weighting  on  T  is  sensible,  since 
if  WDf  W(t,f)  or  some  modified  version  is  to  represent  the  time  behavior  at 
t,  the  transform  on  t  in  line  3  of  (F-6)  ought  not  to  involve  arbitrarily 
distant  values  of T ;  otherwise,  waveform  s  will  then  be  utilized  at  argument 
values  very  different  from  the  time  t  of  interest  and  would  be  non¬ 
representative.  However,  there  is  no  weighting  on  v  in  line  4  of  (F-6), 
thereby  allowing  arbitrarily  distant  argument  values  of  spectrum  S,  from  the 
frequency  f  of  interest,  to  be  considered;  this  unrealistic  feature  of  the 
WDF  is  an  additional  reason  for  its  drawbacks. 

SEPARABLE  SMOOTHING 

If  smoothing  function  V?  is  separable,  then 

V2(t.f)  -  va(t)  Vb(f) 

V2(t,T)  =  Va(t)  Vb(r) 
q2(v,T)  =  Va(v)  vb(T) 

Q2(v,f)  =  Va(v)  vb(f)  . 


Then  (F-l)  gives,  for  example. 
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t  f 

W(t.f)®  va(t)®  Vb(f)  = 

-  y  dT  exp(  -i2irfT)  vb(r)  J  dt1  R(t  -  t'.T)  va(t')  . 

This  has  both  the  desirable  features  of  locally  averaging  the  correlation 
and  suppressing  large-T  contributions .  However,  it  restricts  the  form  of 
averaging  in  the  t.f  plane  and  specifically  disallows  tilted  smoothing 
regions  which  are  not  parallel  to  the  t  or  f  axes. 

GENERAL  GAUSSIAN  TWO-DIMENSIONAL  SMOOTHING 


The  inadequacies  of  smoothing  in  time  alone  or  frequency  alone  suggest 
consideration  of  the  general  two-dimensional  result  in  (F-l): 

A  tf 

W(t.f)*  W(t,f )  ®  v,(t,f)  = 


=  J  dt  exp(-i2irfX)  J  dt'  R(t  -  t'.T)  v2(t’,T)  - 

*  J  dx  exp(-i2trfT)  j'  dt'  s(t  -  t'  4-  |)  s*(t  -  t'  -  |)  v2(t',r)  . 
If  we  let  t1  =  t'  +  J  ,  t2  =  t'  -  2  •  becomes 
W(t,f)  =  JJdt]  dt2  exp(  -i2wf (t-j  -  t2))  s(t  -  t2)  s*(t  -  t,)  v2 


(F-7) 


( F  -8 ) 


F-6 


Now  let  two-dimensional  smoothing  function  have  the  general 
Gaussian  form 

V2(t,f)  =  2/tf  exp{ -a2t2  -  4*2b2f2  -  4wctf)  ,  (F-9) 

where  a,b,c  are  real  constants  and 

2  2  2 

0=  a  b  -  c  .  ( F-10) 


The  scale  factor,  2\Z0\  is  chosen  so  that  the  volume  under  V2  in  the  t,f 
plane  is  1;  this  keeps  the  volume  under  the  smoothed  distribution  In  (F-l) 
or  (F-7)  at  E,  the  energy  of  waveform  s.  In  order  that  V2  tend  to  zero  at 
infinity  in  the  t.f  plane,  we  must  have  Q  >  0.  The  area  in  the  t,f  plane  of 
the  contour  ellipse  of  (F-9),  at  the  1/e  relative  level,  is  (appendix  D) 


Area  =• 


1/2 


1/2 


VaV  -  c2  '  V?' 


( F  — 1 1 ) 


The  transform  on  f  of  Vj  in  (F-9)  is 


v2(t*T)  =  exP 


[-?[ 


2  r2 

Q  1  +  4  +  ictT 


}] 


( F-l 2) 


For  completeness,  the  two  remaining  smoothing  functions  in  (F-l)  are 
q2(v,r)  =  exp jj  ^(aV  +  4*2b2v2  +  4*cf^)]  , 

Q2(v.f)  =  |  exp ^(Of2  + 


We  can  now  determine  the  quantity  necessary  for  evaluation  of  (F-8),  namely 
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*2^“  •  *,  -  l2)  -  ‘F  •»»[-  ^  {<»  *  1  "  (2C>  *?  - 


+  (Q  +  1  -  12c)  t,  «■  2(Q  -  1)  t 


(F  13) 


By  the  discussion  in  (25)  et  seq.,  this  function  in  ( F -1 3 )  is  separable  in 
t1  and  t  if  and  only  if 


Q  -  1 ,  that  is ,  a^b2  -  c2  =  1 . 


(F  -14) 


Then  V2  in  ( F -9 )  Is  a  legal  WDF  and  the  area  in  (F-ll)  becomes  1/2.  Also, 
smoothing  in  (F-9)  is  then  exactly  equivalent  to  the  Gaussian  smoothing 
considered  previously  in  (147). 


We  are  Interested  here,  however.  In  the  more  general  case  of  V?  where 
Q  is  not  necessarily  1,  and  therefore  V?  is  not  a  legal  WDF.  If  we 
substitute  ( F-l 3)  in  (F-8),  the  smoothed  WDF  becomes 


W(t,f) 


x(  t-j  )  x 


*(t2)  exp 


Oil  t  t 
2b^  1  1 


(F  IS) 


where 


r  t2 1 

xft^  =  s  (t  -  tj)  exp  I  ■ -i 2irf t1  -  (0  *■  1  +  i2c)  — 2J  . 
By  expanding  the  exp  in  (F-15)  in  a  power  series,  there  follows 

a(t.f)  -  Fl  ^  (~?)"  |[‘lt1  «<v  *?r  • 


(F  -16) 


(F-l 7) 


1R  8 22b 


L't^lL'1.1  .h' Ji* Ji».ha J|'  k* lii  L« 


I**  >v  ttftU'lKtiiMt.- 


It  is  obvious  from  ( F  -17)  that  a  sufficient  condition  for  smoothed  WDT  W 


to  be  non -negative  is  Q  <  1,  that  is 


i^b^  -  c^  <  l 


( F-l 8) 


(lhe  special  case  of  c  -  0  was  given  in  [11.  (5)].)  When  this  condition  is 
used  in  ( F -11 ) ,  we  see  that  the  area  of  the  concentration  ellipse  of  (F-9)  is 


Area  >  1/2 


(F-19) 


Thus,  smoothing  with  the  Gaussian  two-dimensional  function  V?  in  (F-9)  always 

results  in  a  non-negative  distribution,  provided  that  the  area  of  the  ellipse 

at  the  1/e  relative  level  is  greater  than  or  equal  to  1/2.  It  is  not 

necessary  that  V?  in  (F-9)  be  a  legal  WOT;  that  is,  the  area  of  the  ellipse 

need  not  be  precisely  1/2.  However,  the  most  concentrated  V,,  in  (F-9),  that 

guarantees  nonnegative  results,  has  area  1/2.  The  only  restrictions  on 

2  2  2 

parameters  a.b.c  are  given  by  0  <  0  <  1 ,  that  is,  0  <  a  b  -  c  <1. 


F-9/F-10 
Reverse  81ank 


I 
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APPENDIX  G.  DERIVATION  OF  MINIMUM-SPREAD  WDF 

The  short-term  spectral  estimate  of  waveform  s,  relative  to  weighting  u, 
is  given  by  the  double  convolution  (106), 


tf 

ws(t,f)  $  wu(t,f)  , 


(G-D 


of  the  WUFs  of  s  and  u.  It  is  therefore  important  to  use,  for  weighting  u, 
a  function  which  has  as  narrow  a  WDF  as  possible,  so  that  the  smearing 
implied  by  (G-l)  is  minimized.  In  particular,  since  we  are  interested  in 
analyzing  waveforms  with  linear  frequency  modulation,  we  are  interested  in 
inimizing  the  spread  of  WDF  Wu>  as  measured  by  the  quadratic  quantity 


mi 


I  =  JJdt  df  Wu(t,f)  (f  -  0ct)2  , 


where  0c  is  a  specified  (observed)  slope  in  the  t.f  plane. 


By  expanding  the  quadratic  in  (G-2),  we  obtain  spread 


where 


1  “  lo  *  !1  *  *2  • 


l0  -  JJdt  df  f2  Wu(t,f)  , 

I)  =  -26c  ff  dt  df  1  f  Wu(l’f>  • 

\2  -  B2  dt  df  t2  Wu(t,f)  . 


(G-2) 


(G-3) 


(G-4) 

G-l 


Reference  to  (43),  (28),  and  (12),  respectively,  allows  the  terms  in  (G-4) 
to  be  simplified  and  expressed  solely  in  the  time  domain  as 

i  r  2 

^  "  2  Sdt  lu'(t)i  • 

0  4/  J 

^  ^  |dt  t;  Imfu'(t)  u*(t)}  , 

l2  =  0*  f  dt  t2!u^)l  •  (6 -5) 

Adding  these  results  together,  the  spread  in  (G-3)  becomes 

I  «  fdt  |u'(t)  -  iact  u(t)|  2  .  (G-6) 

4i 

where  we  define 

a.  =  2*Gc  •  (G-7) 

Observe  that  the  spread  I  in  (G-6)  is  nonnegative,  for  all  weightings  u. 


The  function  that  minimizes  I  in  (G-6)  is 


u(t)  =  aQ 


aQ  complex, 


(G  8) 


for  which 


J...V  .»t-VAV,-  JSJA, t£h!  .tJ-Mrm 


^<4  W  »  V  I  .->1^ 
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I  -  0  ,  Wu(t,f)  =  |aop  S(  f  -  Bct) 


(G-9) 


That  is,  the  WOP  is  concentrated  on  the  f  =  8  t  line  in  the  t,f  plane. 

c 

However,  the  energy  of  (G-8)  is 


Eu  =  06  . 


( G— 1 0 ) 


which  is  unacceptable. 


If  we  attempt  to  approximate  (G-8)  by  unit  energy  weighting 


?  "1/4  /  t2  J 

u(t)  =  (*</)  exp/-  +  i  f  XT k 


( G— 1 1 ) 


the  spread  of  u,  as  given  by  (G-6),  turns  out  to  be 


i  1  r  1  .  2/  v2l 

'  ‘  e7 [j  °  (“  ’  •«'  J  • 


(6-12) 


Now  if  a  =  <*c,  that  is,  the  linear  frequency  modulation  parameter  a  in 
weighting  u  is  exactly  equal  to  given  quantity  <*c  (from  (G-7)  and  (G-2)), 


the  spread  is 


*  “  2  2  >  ®  ' 
8ir  <j 


(6-13) 


However,  as  duration  o  of  weighting  (G-ll)  gets  larger,  the  spread  I  tends 
to  zero,  even  though  the  weighting  has  finite  (unit)  energy.  Also,  then 
(G-ll)  tends  to  a  scaled  version  of  (G  8)  at  each  fixed  t. 


ft 

1 


m 

% 


K 

k> 


I 

% 

a? 

Lei 
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In  order  to  eliminate  these  undesirable  features  of  the  weighting,  it 
will  not  be  sufficient  to  minimize  spread  1,  subject  only  to  a  constraint  on 
the  energy  of  u.  Rather,  it  will  also  be  necessary  to  constrict  the  time 
duration  of  the  weighting  u.  Accordingly,  we  will  minimize  spread  I  in 
(6-2) -(G-6),  under  the  two  constraints  that 


fdt  ju(t)|  2  =  1  =  JJdt  df  Wu(t,f)  , 


Jdt  t2  |u(t)|  2  =  ~f  =  f J 


dt  df  tc  Wu(t,f)  ; 


( G— 1 4 ) 


see  (M)  a..i  (12). 


Thus  consider 


0  =  J  dt  |u'(t)  -  1act  u(t)|  2  +  X  Jdt  |u(t)|  2  ♦  u  |  dt  t2  |u(t)|  2  ,  (G-15) 


where  X  and  v  are  real  Lagrange  multipliers.  Replacing  u(t)  by 
u(t)  +  en(t),  where  n(t)  is  an  allowed  variation,  we  have 


0  4-  *Q  =  Jdt  [u'(t)  4-  cn'(t)  -  1act  (u(t)  4-  cn(t) )  ]  [u  '*(  t)  +  c%'*(t)  4- 


4-  iact  (u*(t)  4-  c%*(t))]  4-  Jdt  (u(t)  4-  cn(t)  ]  [  U*(  t )  4-  cV(t)]  (X  4-  J1  t?) 


(G-16) 


G-4 


rw*r*vT« "jiM"jrKjr*3r*jFjinji-*  -jiujip^’vtw  wi#rn* rwnww r*w*  wwf- 
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Since  the  coefficient  of  e*  must  be  zero,  In  order  for  u(t)  to  be  the 
optimum  [15],  we  require  the  following: 

(dt  n'*(t)  [u ' (t)  -  iact  u(t)]  4-  j*dt  n*(t)  [u'(t)  -  iact  u(t)]  iact  + 

4-  Jdt  n*(t)  u(t)  (X  4-  yt2)  -  0  for  all  n(t)  .  (G-17) 

We  now  integrate  by  parts  on  the  first  term  of  (G-17),  and  presume  that 
n.u.u'  all  decay  to  zero  at  ±tt>.  Since  n  is  arbitrary,  its  coefficient 
under  the  integral  must  be  zero;  namely,  we  find  that  u  must  satisfy  the 
following  differential  equation: 

-  u  (t)  4-  i2act  u * ( t )  4-  (iac  +  a2t2  4-  \  +  yt2)  u(t)  =  0  for  all  t  .  (G-18) 

If  we  try  solution 

u(t)  -  a  exp^  ct^  •  a.c  complex  ,  (G-19) 

in  (G-18),  we  find  that 

t2(-c2  4-  12acc  4-  o2  4-  y)  -  c  4-  i«c  f  X  *  0  for  all  t  .  (G-20) 

Ihen  a  solution  of  the  form  (G-19)  exists  with  the  choices 

c  -  iac  t  f?  ,  x  =  t  /T.  (G-21 ) 


6-5 
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lo  determine  u(t)  explicitly,  we  substitute  (G-21)  into  (G-19)  to  get 

u(t)  -  a  exp^i  ^  ac  t2  ±  t2^  .  (G-22) 

When  the  two  constraints  in  (G-14)  are  satisfied,  there  follows,  for  the 
optimum  weighting. 


u(t) 


('  «t)  exp  (  ^2  *  1  "2  t?)  • 


( G— 23 ) 


This  is  linear  frequency  modulation  with  a  Gaussian  envelope. 

The  minimum  value  of  the  spread  1  in  (G-6)  for  the  optimum  weighting 
(G-23)  is 

1 


minimum  1  = 


a  2  2  ’ 
Sir  oc 


(G-24) 


and  the  corresponding  WOF  is 


Wu(t,f)  -  2  expj  - 


=  2  exp  - 


[  ^  *  o  *  4irft®c  -  ■ 


( G-25) 


where  e  =  a  o  .  The  area  of  the  contour  ellipse  at  the  1/e  relative 
c  c  c 

level  is  1/2  in  the  t,f  plane. 


Ihe  mean-square  time  extent  of  the  optimum  weighting  u  in  (G-23)  is 

p 

o‘/2,  as  required  by  constraint  (G-14).  The  mean-square  frequency 
extent  is  obtained  from  voltage  spectrum 


1  +  e 


o  2  2 

8*  a 

c 


(G-27) 


G-7/G-8 
Reverse  Blank 


APPENDIX  H.  EXAMPLE  OF  WDF  PROCESSOR 

The  processor  of  interest  here  is  depicted  in  figure  7,  while 
representative  characteristics  of  the  waveforms  and  devices  are  sketched  in 
figure  8.  The  mean  output  is  given  by  (173)  and  the  variance  is  given  by 
(188).  We  will  use  the  definitions  and  results  in  (165)  (189)  freely  in  the 
fol lowing. 


INPUT  1NF0RMA1 ION 


The  Input  signal  waveform  to  figure  7  is  Gaussian -modulated  linear 
frequency  modulation: 


s0(t)  *  a. 


o  ♦ 1  i  t?J  •  kl2«o  • 


(H-l) 


where  aQ  can  be  complex.  The  instantaneous  input  signal  power, 


k(t)l2  -  kl2  «*p<-  t2/«o>  • 


(H-2) 


peaks  at  t  =  0  and  has  effective  duration  a  . 

o 


The  corresponding  signal  voltage  density  spectrum  is 


7R  8226 


The  energy  density  spectrum  is 


2*Ja  I 2 a2 


.22  ,2 
4*  a  f 
0 


(H-4) 


fm 


cv 

I 


which  peaks  at  f  =  0.  When 


f  =  ±  f  =  i  - 


1  f  0" 

_ g 

2*o 


(H-5) 


the  energy  density  spectrum  Is  reduced  to  1/e  of  its  peak  value;  hence  f 
Is  a  measure  of  the  bandwidth  of  the  linear  frequency  modulation  waveform. 


The  input  noise  nQ  to  figure  7  has  a  white  spectrum 


Gn  (f)  =  Nd  for  a11  f  • 
o 


(H-b) 


£ 


\* 

> 


The  filter  transfer  function  is 


H(f)  -  exp 


(3- 


(H-7) 


which  peaks  at  f  =  0;  this  coincides  with  the  signal  spectrum  peak, which 
means  that  we  are  considering  the  most  fortuitous  situation.  The  weighting 
in  figure  7  is  taken  to  be 


v(t)  -  exp  I 


(H-8) 


vr  c~  n»  kji  mmitTL*  i.w  v  Tvrrrv  rows  v*.-rw.  vw.-giuw  ■wi*rj»rj«  rm  w 
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which  peaks  at  t  =  0;  again,  this  coincides  with  the  signal  peak  and  is  most 
favorable.  The  maximum  values  of  H  and  v,  being  equal  to  1,  are  chosen  for 
convenience,  without  loss  of  generality;  absolute  level  does  not  influence 
t fie  performance  of  the  processor  in  figure  7. 

CALCULATIONS  OF  BASIC  FUNCTIONS 


We  will  make  frequent  use  of  ( 68 ) -(70)  in  evaluating  the  following 
quantities  which  are  needed  in  (165)  -(189);  the  choices  for  Gaussian 
functions  for  s0,H,v,  above,  were  made  for  analytic  simplicity,  since  the 
various  integrals  can  be  conducted  In  closed  form.  More  general  cases  would 
require  numerical  integrations. 

The  noise  spectrum  at  the  output  of  the  filter  is 

Gn(f)  -  Gn  ( f )  |h ( f  )|  *  Nd  exp(  -  f  2/B2)  .  (H-9) 

o 

Ihe  corresponding  noise  correlation  is 

Cn(x)  =-  |  df  exp(i2wfx)  Gn(f)  *  H  Nd  B  exp^wVx2)  .  (H-10) 

The  auxiliary  spectrum  in  (178)  is 

6n?)(f)  '  f  dr  exP(  -Hirfr)  C2(x)  =^|’ N2  B  exp(-  2f2/B2)  .  (H-ll) 
The  instantaneous  correlation  of  the  weighting  is 


H-3 


» 


3 


i 


(H-12) 
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:,I)  -  v(t  ♦  |)  v*(t  -  f)  -  exp ^  ^  -  -~y 


and  its  corresponding  WOf-  is 


Wv(t,f )  =  |  dT  exp(  -i 2*fx )  Rv(t.t)  =  2  y?  L  exp  ^  —  -  4ir2L2f^  .  (H-13) 


The  filter  impulse  response  is 


h(T)  =  Jdf  exp( i2*fT)  H(f)  -  B  exp(  -2*2B2r2)  , 


leading  to  filter  output  signal 


s ( t )  -  [  dr  h(r)  sQ(t  -  t)  -- 


4^),/2  4  ^2r^V2]- 


where 


0A  -  <*  aA  ,  0  -  (2»Ba 

0  0  0  0 


(H  -14) 


(H-l  5) 


(H  16) 


This  filter  output  signal  is  again  a  Gaussian -modulated  linear  frequency 
modulation  waveform. 


In  general ,  for  signal 

s ( t )  -  cq  exp(  -ct2/2)  ,  cQ,c  complex  ,  E  =  WfcQ|  2//cT7f  (M-l 7) 


v.  9* ♦  v  •  .v  *  rryT!v*7 
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the  WOF  is 


W.(t.f)  -  2 fir-=—  exp  - 

K  L 


t2  jc i  2  4-  4w2f  2  4-  4*f  tc  • 


(H-18) 


When  applied  to  example  (H-15),  we  identify 


f  n  V/2 
o  ao\l  4-  o  -  ieoj  * 


,  2n2  1  10o 

c  =  4* 8  rrr^ 


(H  -19) 


to  obtain 


W2-  W2*.  M2-«**>4 


o  o  di  0  0  to 

„  ,  2d2  1  „  .  2d2  _ 0 

Cr  r  4*  8  *  Ci  ‘  -  4*  B  T2 


1  4-  e 


(H-20) 


where 


D1  =  1  4-  D  4-  e2  ,  02  -  (1  4-  D)2  f  e2 


( H— 21 ) 


When  substituted  in  (H-18),  there  follows,  for  the  WOF  of  s  in  (H-15), 


ws(t,f) 


>(sf)  «»[-  5;["  -  *  °i(g)  -  "..«««(£ 


(H-22) 


As  bandwidth  8  of  filter  H  in  (H-7)  tends  to  infinity,  then  D  -» <*, 

01  -»  0,  02  -»  02,  and  (H-2?)  yields 

W$(t,f)  -4  2E0  exp (1  4  e2)  ^2  -  (2woQf )2  4-  4*eQftjas  B  -»<*  .  (H-23) 


*  k’.'SV.V 
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This  agrees  with  (91)  and  (78).  More  generally,  in  order  to  keep  the  scale 
factor  (D/0^)1^2  in  (H-22)  near  1,  we  need 


0  >  1 


that  is,  B  > 


fH 

2"o  ' 


(H-24) 


according  to  ( H -2 1 )  and  (H-16).  But  this  latter  quantity  is  just  the 
bandwidth  of  signal  sq;  see  (H-5).  Thus  condition  (H-24)  yields  the 
physically  intuitive  statement  that  the  filter  passband  should  be  wider  than 
the  input  signal  bandwidth,  in  order  not  to  decrease  the  peak  value  of  the 
WDF  of  the  filter  output  signal. 

The  area  of  the  elliptical  contour  of  the  general  WDF  in  (H-22),  at  the 

1/e  relative  level,  is  1/2  in  the  t.f  plane,  regardless  of  the  values  of  any 

of  the  parameters  of  the  input  signal  and  filter;  this  follows  by  the  direct 

use  of  (0-1),  (0-19),  and  (0-20).  It  is  also  consistent  with  the  general 

fact  that  this  is  true  for  any  signal  of  the  form  of  (H-17),  as  may  be  seen 

by  application  of  appendix  0  directly  to  (H-18),  where  c  and  c  are 

o 

arbitrary  complex  constants. 

The  peak  height  of  the  signal  WDF  in  (H-22)  is  2Eo^D7iTJ;  hence,  the 

product  of  peak  height  and  effective  area  is  e0V57E7,  wh  ich  is  just  the 

energy  of  s: 

E  =  |  dt  |$(t)|  2  =  jjdt  df  Ws(t,f)  .  ( H  -25) 


This  follows  directly  from  (H-15)  or  (H-17). 
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ihe  tilt  of  the  major  axis  of  the  elliptical  contours  of  (H-22)  is  given  ; 

by  0  radians  in  the  (2«Bt,f/B)  plane,  where  i 

i 

e  ! 

tan(20)  *  .  (H-26) 


according  to  ( D -5 )  and  (H-21).  ■ 

* 

I 


MEAN  SIGNAL  0U7PU1 


The  mean  signal  output  of  the  WDF  processor  in  figure  7  is  given  by 
(171)  as 


l 

t 

i 


I 


a  -  J dr  exp(  -i2«fx)  Ry(t,T)  Rs(tfr)  = 

=  jdt  exp(-i2wfT)  R~(t,T)  =  W?(t,f)  ,  (H-27) 

where 


s(t)  --  v(t)  s ( t )  «  cQ  exp( -c  t2/2) 


( H -28 ) 


and 


c  =  c  + 


(H-29) 


according  to  (H  8)  and  (H -1 5) -(H-l 9) .  By  analogy  to  (H-18),  we  have 


H-7 


I 


,1'tJS  i-%  It.  lWftVlU  lU'tVl 
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Jc  I  2  ["  t2  |  C  J 2  f  4-»r2  f 2  f  4irf  t  c. 

a  -  W  (t.f)  --  2  V?  exp  -  - = - 

?  tfc  Cr 

'  r  L_ 


nQm  exp  ^[h2  X2  f  0?y2  -  2RD0 
\  1  /  L.  1 


o  xyjj  • 


( H -30) 


where  we  def i ned 


R  -  ( 2irBl;  ,  x  -  t/L  ,  y  =  2«fL  , 


=  02  f  R0]  =  (1  +  0) ( 1  <-  D  +  R)  4-  (1  *  R)  , 


H  =  D0  4-  2R0  f  R2(l  +  02)  -  (1  ¥  0  4-  R)2f  (1  4-  R)2©2  ,  (H-31) 

c.  c  I  C  0 


and  used 


H  R0©  „  H„ 

~  _  1  ~  _ 0  1 2  2 

cr  "  2  *  Ci  =  2  •  !ci  4 

L  02  rO?  l  02 


( H-32) 


The  area  of  the  concentration  ellipse  of  W~  in  (H-30)  at  the  1/e 


relative  level  is  1/2  in  the  t.f  plane,  regardless  of  the  sizes  of  D  and  R; 


so  the  signal  WOf  is  not  spread  by  the  filtering  and  weighting  operations  in 


figure  7,  at  least  for  this  example  of  (H-l)  coupled  with  (H-7)  and  ( H -8 ) 


The  peak  height  of  WOF  W^  is  given  by  the  leading  factor  in  (H-30) 


Since  the  effective  area  of  this  WDF  is  1/2,  the  product  of  peak  height  and 


effective  area  is  E  ^DR/H^,  which  is  just  the  energy  of  s: 


t  .  Jdt  |?(t)|  ?  =  i  f - 26 - Y 

-  °V<>  *  0)<1  +-D+R)  +-(1  >  R)  8  ‘J 


( H  -33) 
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Ihe  parameters  eo,D,R  are  given  in  ( H  -1 6 )  and  (H-31). 


We  now  define 


1  +  R 


L  1  1  .  r>>  2 

-  •••  =**  ;  —  =  —  ¥  ( 2irB) 

/l  +  (2wBL)  L 


( H  —34 ) 


Then  (H-30)  and  (H-31)  become 


¥*•'>  '  ^  *  02(HR)(2,Lbf)2  -  2RDV*V>(l$] 


( H -35 ) 


Ihe  major  axis  of  the  elliptical  contours  is  at  angle  If  radians  in  the 


(t/l_b,2irLbf )  plane,  where 


0(1  +  R) 

tan(2^)  *  1  +  R  +  0  +  RD/2  • 


( H —36 ) 


Given  measurements  or  observations  B  and  oq,  this  can  be  immediately 


solved  for  eQ,  where  it  is  presumed  that  B  and  L  are  known  since  they  are 


under  our  control . 


As  alternative  fundamental  parameter  is  more  useful  than  the  above;  we 


introduce 


M  =  (L/o0)2  , 


( H  —37 ) 


which  is  the  square  of  the  ratio  of  weighting  time  to  input  signal  duration. 


Then  (H-31)  yields 


R  -  DM  , 


(H-38) 


a 


a 
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where  0  -  (2irBoQ)  just  as  defined  in  (H-16).  Eliminating  R  in  favor  of 


0  and  M,  the  peak  height  of  W„  in  (H-30)  becomes 


(1  +  D)(l  fO  ^  OH)  t  (1  r  DM)0o 


As  checks  on  this  quantity,  observe  that  as  D  ->  *> ,  the  factor  of  PE 

1  /2 

in  (H-39)  approaches  (M/(M  +  1))  ;  in  order  to  keep  this  latter  quantity 

near  1,  we  need  M  >  1,  that  is,  L  >  o  .  This  is  consistent  with  physical 

reasoning  on  figure  7.  Alternatively,  as  M  -»o t,  the  factor  of  2Eq  in 

2  1  /2 

(H-39)  approaches  (0/(1  +-  0  +■  eQ))  ;  in  order  to  keep  this  near  1,  we 
need  (H-24)  to  be  observed,  just  as  before. 


More  generally,  in  order  to  keep  the  signal  factor  in  (H-39)  near  unity, 
we  need  to  choose  the  combination  of  D  and  M  large  enough.  For  example,  to 
keep  the  factor  at  value  F,  that  is,  maintain 


=  F  (<  1)  , 


(H-4C 


(1  +  D)(l  +  0  +  0M)  +  (1  -t-  OM)0 


we  need  to  choose  L  in  (H-37)  such  that 


F2[(l  i-  D)2  *  e2] 
D[D  -  F2( 1  f  0  f  e2)] 


However,  this  relationship  is  useful  only  for 


D  >  — £ — ~(1  +  e2)  . 
1  - 


(H  41 


(H-42 
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A  representative  sketch  of  (H-40)  is  displayed  in  figure  H-l .  Small  values 
of  parameters  M  or  0  are  not  realizable  without  the  other  parameter  tending 
to  infinity,  in  order  to  maintain  the  factor  in  (H-40)  at  value  F.  Three 
numerical  examples,  for  9q  =  0,1,5  are  given  in  figures  H-2,  H-3,  and  H-4, 
respectively.  The  larger  F  values  can  only  be  achieved  through  rather  large 
D  and/or  M  values. 


lhere  is,  however,  a  minimum  value  of  the  product,  MD, 
realize  a  specified  value  F  for  the  factor  in  (H-40),  when 
specified.  In  fact,  we  find  from  (H-41)  that,  for  given  F 
product  MO  is  minimized  by  the  choices 


required  to 
9o  is 

and  0  ,  the 
o 


opt 


o  J 


( H— 43) 


The  value  M  .  is  relatively  insensitive  to  e  ;  in  fact,  it  varies  from 
opt  J  o’ 

4F2/(1  -  F4)  to  2F 2 ( 1  F2)  as  eQ  varies  from  0  toad,  which  is 

less  than  a  2:1  variation. 


The  corresponding  minimum  product  is 

.2 


(MO) 


min 


(i  -  f2) 


2.2 


iT 


♦  f  tV(it  e:)()  ♦  f 


^  • 


(H-44) 


H-l  1 


2 


4 


lo 


12 


)4 


Figure  H-3.  Plot  for  0O  =  1 


32  4?  p  £4  ro 


Figure  H-4.  Plot  for  0O  =  5 
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H 


I 


In  particular, 


‘""’min 


0  -  f2) 


-  (l  4-  f  eo)  4-  l  as  e0  -» ot>  ; 


( H  -45) 


in  fact,  this  is  a  good  approximation  except  near  eQ  -  0.  Thus,  large 
amounts  of  linear  frequency  modulation,  or  values  of  F  near  1,  require  very 
large  MD. 


At  the  other  extreme, 


for  e  =  0  . 
2  o 


( H  -46) 


For  example,  if  F  =  1/V?,  this  product  is  8;  thus  relatively  large  values  of 
the  product  are  required,  even  at  the  low  end  where  there  is  no  linear 


1R  B?25 


As  a  particular  numerical  example,  for  0Q  =  0,  F  =  1/| fT,  we  find 


°opt  =  3  *  "opt  =  3  '  (MD)min  =  8  ' 


( H  -4  7 ) 


So  both  0^^  and  MQp^  are  somewhat  larger  than  unity,  even  for  F  equal  tt> 
1/V?.  All  these  conclusions  are  drawn  relative  to  the  mean  signal  output 
alone;  we  now  consider  the  noise  output  contributions. 


MEAN  NOISE  0UTPU1 


The  mean  noise  output  of  the  WDF  processor  in  figure  7  is  given  by 


(172),  (H-9) ,  (FI-13) ,  and  (H-31)  as 


b  =  Wy(t.f)  «  Gn(f)  = 


.  N  /_jlV/2  expr  t£  _  d  _r_1 

Nd  (l  +  r)  6  P  l2  g2  1  +  R  J  ' 


( H  -48) 


(As  L  -»c0,  that  is,  no  weighting,  then  R  -»«*,  and  b  -»  N^  exp( -f  /d  ) 
as  expected.) 


Gn(f). 


The  noise  factor  in  (FI-48),  namely 


/_R_  V/2  _  /  DM  V 
+  r /  vi  +  dm; 


(FI -49) 


is  virtually  unity  when  the  mean  signal  degradation  is  small,  according  to 
the  results  of  figure  H-5.  Thus  the  ratio  of  peak  signal  to-noise  means  is 
approximately  2Eq/N^,  according  to  (H-30)  and  (FI-48).  These  latter 


»  V'  ,  V,  V  |*»  V  :  V*  V*1  I  V  *1*  1  *  *  I  4*»  I  ’?  i 


*v  V  y  *  V\ 


. 1  t  «  •  v 


quantities,  Eq  and  N^ ,  are  directly  the  input  parameters  to  the  WDF 
processor  in  figure  7;  see  (H-l)  and  (H-6). 

VARIANCE  OF  NxN  TERM 


The  variance  of  the  NxN  term  at  the  WDF  processor  output  is  given  by  the 
third  line  of  (188)  as 


VNN  =  2  W2(t,f)  0  &<2)(f) 


R 

2t2  2f 2  R 

vr - <  e*P 

yl  4-  R 

L  L2  '  B2  W  RJ 

(H  50) 


Here,  we  also  used  (H-ll)  and  ( H -1 3 ) .  As  l  -♦  <*,  then  H  -»<*>,  R  -»«>,  and 
V.,..  ><*>.  Alternatively,  as  B  -*<*>,  then  0  -♦  <*>.  R  -*<»,  and  -»«*>. 

NN  Hfi 

These  results  for  this  particular  example  confirm  the  general  observations 
in  the  sequel  to  (188) . 


Ihe  standard  deviation  of  the  NxN  term,  namely  \/v  ',  is  precisely  equal 

to  the  noise  mean  output  in  (H-48)  for  all  t,f,  except  for  a  constant  factor 
(1  4-  R)1/4  -  (1  4-  OM)1^4.  Also,  the  axes  of  the  elliptical  contours  of 
(H-48)  and  (H-50)  are  parallel  to  the  t  and  f  axes  and  are  independent  of 
eQ,  the  amount  of  linear  frequency  modulation  in  the  input  signal. 
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VARIANCE.  OE  SxN  URM 


The  variance  of  the  SxN  terms  at  the  WOE  processor  output  is  given  by 
the  last  line  of  (188)  or  by  double  the  results  in  (182).  Upon  substitution 
of  ( H  -1 0) ,  (H-12),  and  (H-15)  in  (182)  and  an  extreme  amount  of 
manipulations,  there  follows  variance 


4E  N. 
o  d 


exp[-£(t,f ) ]  . 


( H -5 1 ) 


where 

H-  =  ( 1  ♦  D  R/2) ( 1  +  0  ♦  R  f  OR/2)  f  (1  ♦-  R)(l  *  R/2)  02  (H-52) 

and  £(t,f)  is  an  elliptical  function  with  minimum  value  at  t  =  f  =0. 

Namely, 

,t2  o  *  2  .  *2  »  ,  ■  2 

f(t,f)  .  21^  «.  ^t2  -  .  (h-s3> 

L2  2(|a|Z-YZ) 

where 


c 


2  2 

4*  8 


ie_ 


l  *  D  -  ie_ 


«  *  x(c  *-  4  f  2  *2B2)  , 

L 

u  -  -  (f  t  f  i2*f )  , 


1  2„2 


/  ii  r  *  v 


The  quantity  in  the  denominator  of  ( H  S3 )  can  be  simplified  to 


2<H2  -  T2)  =  -r-  ;  ("-Ml 

2L,D? 

however,  £(t,f)  has  not  been  reduced  lo  its  most  compact  form, 

b}t2  <-  b2f2  f  2b3tf  ,  (H-b6) 

due  to  the  excessive  amount  of  labor  required  to  simplify  and  obtain 
bl ,b2’b3' 

QUALITY  MEASURE  OF  PERFORMANCE 


We  define  a  quality  measure  for  the  WOF  processor  output  in  figure  7  as 


Difference  of  mean  outputs  _ 
Standard  deviation  of  output  ~ 


SN 


f  V 


nn) 


w 


(H-S7) 


The  relevant  quantities  are  given  by  (171)  and  (188)  generally.  For  the 
specific  example  in  this  appendix,  the  quality  measure,  at  peak  signal 
location  t  =  f  =  0,  is  obtained  by  combining  (H-30) , (H-50) ,  and  ( H -5 1 ) : 


ovr »  r' 

JUi _ 


( H-58) 


For  convenience,  we  repeat  the  parameter  definitions  here: 
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®o  :‘.V0'  <2*B°0)2  •  M  *  iu”o'>2  ■ 

R  =  DM  . 

H]  *  (1  *  D)  ( 1  «•  D  R)  «■  (1  *  R)®o  . 

H3  =■-  (1  ♦  D  f  R/2) ( 1  +  D  +  R  *  DR/2)  +  (1  +  R) ( 1  +  R/2)e*  .  (H-59) 

It  has  already  been  observed  in  (H-24)  and  in  the  sequel  to  (H-39)  that 

D  >  1  f  0Q  and  M  >  1  are  desirable,  in  so  far  as  the  mean  signal  output  is 

concerned.  However,  if  filter  bandwidth  B  (D)  is  made  too  large,  then  too 

much  noise  is  allowed  through;  alternatively,  if  weighting  duration  L  (M)  is 

made  too  large,  a  noise  degradation  also  results.  Thus,  it  is  expected  that 

2 

the  quality  ratio  0  will  peak  for  D  in  the  neighborhood  of  1  +  0Q  and 
for  M  near  1 . 

It  should  be  observed  from  (H-58)  that  even  if  input  signal-to-noise 
measure  Eq/N^  gets  extremely  large,  the  quality  measure  0  behaves 
according  to  \jl  /N  ’  and  not  E  /N..  This  is  due  to  the  saturation 
effects  caused  by  the  SxN  terra  in  the  denominator  of  definition  (H-57);  it 
can  also  be  seen  directly  from  the  quantitative  result  In  (H-51),  where 
variance  Vcu  is  directly  proportional  to  input  signal  energy  E  as  well 
as  the  noise  density  level  N^. 


The  quality  ratio  Q  in  (H-58)  is  plotted  versus  H  in  figures  H-6,  HI, 


H-8  for  9  -  0,1,5,  respectively.  The  input  ratio  t  /N,  is  kept  at 

o  o  0 

value  20  in  all  cases;  the  only  other  fundamental  parameter,  0,  is  varied 

over  a  range  wide  enough  to  encompass  the  maximum  of  Q.  However,  for  ease 

of  plotting  the  results,  the  values  of  0  which  are  less  than  the  critical 

value,  which  leads  to  the  peak  Q,  are  separated  from  those  that  are  greater 

than  the  critical  value.  For  example,  in  figure  H-6,  D  =  1  leads  to  the 

maximum  value  of  Q  that  can  be  achieved  for  any  value  of  M;  thus,  the  upper 

part  of  figure  H  6  contains  results  for  D  <  1,  while  the  lower  part  contains 

the  remainder  for  0  >  1.  Ihe  corresponding  critical  values  of  0  are  8  and 

80  in  figures  H-7  and  H-8,  when  e  =  1  and  5,  respectively. 

o 

One  important  observation  that  is  made  apparent  by  these  figures  is  that 

near  the  maximum,  the  quality  ratio  0  is  not  too  sensitive  to  M  and  D;  that 

is,  the  maximum  is  broad  in  the  neighborhood  of  the  best  parameter  pair 

M,0.  It  should  also  be  observed  that  as  0Q  increases,  the  peak  value  of  Q 

decreases,  although  the  decrease  is  not  very  significant,  at  least  over  the 

range  9q  -  0,1,5  used  here.  Finally,  the  values  of  the  peaks  in  these 

figures  are  slightly  less  than  \JT~7n7,  as  anticipated  above. 

*  o  d 
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APPENDIX  I.  SMOOTHED  WDF  FOR  s(t)  =  t  exp(-t  /2) 


For  the  waveform 


s(t)  =  t  exp(-t  /2)  for  all  t  , 


(1-D 


the  WDF  is 


Ws(t,f)  =  2 fir'  exp( -t2  -  4*2f2)(t2  +  4ir2f 2  -  ])  =  2/7  exp(-r2)(r2  -  j)  ,  (1-2) 


with  energy 


=  Jdt|s(t)l 


=  fiT/2 


d-3) 


Contour  plots  of  the  WDF  in  (1-2)  are  concentric  circles  in  the  (t,2irf) 
plane;  in  fact,  (1-2)  is  a  function  only  of  r2  =  t2  +  (2irf)2.  The 


origin  value  of  W^  Is  -2E  =  -/7,  and  the  WDF  is  negative  for  r  <  1//?, 
while  it  is  positive  for  r  >  !/)[?. 


Let  us  smooth  this  WDF  with  the  most  compact  WDF;  namely,  use  the 
Gaussian  weighting  function  in  (G-23)  with  WDF  (G-25)  with  ec  r  0,  =  1 


Wu(t,f)  -  2  exp(-t2  -  4,r2f2)  =  2  exp(-r2)  .  (1-4) 


The  reason  for  these  parameter  choices  of  ©c  and  o is  that  the  contours 
of  (1-4)  are  also  circles  in  the  (t,2wf)  plane  and  exactly  match  those  of 
the  waveform  WDF  in  (1-2);  this  should  lead  to  minimal  spreading. 


b 
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The  result  of  smoothing  (1-2)  by  (1-4)  is 

112  tf 
|su(t,f)p  =  ws(t,f)  S>  wu(t,f)  - 

1,—  ..  2  .  2,2,  /  1  . .  2  .  2,2,\  lj-,  2  .  2.,,  /T  ,, 

*  4>i^(t  +  4ir  f  )  exp^-  2<t  +  4ir  f  y  =  r  exp(-r  /2)  ,  (1-5) 

which  has  volume  E  as  given  by  (1-3).  Again,  this  is  a  function  only  of 
2 

r  ,  but  it  is  never  negative.  This  smoothed  distribution  is  zero  at  r  =  0, 
and  peaks  at  r  =  V?  with  value  .326.  By  contrast,  the  WDF  in  (1-2)  is  -1.77 
at  the  origin,  a  large  negative  value.  However,  the  waveform  WDF  in  (1-2) 
does  decay  faster  than  the  short-term  spectral  estimate  in  (1-5);  this  is  an 
example  of  the  tradeoffs  that  must  be  accepted  when  using  short-term 
spectral  estimation  versus  the  WDF. 

To  lend  credence  to  (1-5)  as  a  better  measure  of  the  time-frequency 
content  of  s(t),  we  observe  that  at  t  =  0,  the  center  of  gravity  of  (1-5)  is 


Then  we  expect  that 

A  sin(2wfQt)  (1-7) 

ought  to  be  a  good  fit  to  s(t)  of  ( I -1 )  for  t  near  zero.  In  fact,  plots  of 
( I  -1 )  and  (1-7)  for  A  =  exp(  — . 5 )  overlap  for  -1  <  t  <  1. 


1-2 


If  we  attempt  this  same  procedure  for  WDF  W$  in  (1-2),  the  denominator 


is 


o o 

f  df  W$(0, 
Jo 


f)  -  0  . 


(1-8) 


giving  rise  to  f  *  • ,  which  is  useless. 


With  respect  to  t  =  1  instead,  we  find  center  of  gravity 


¥> 

Xdf|su(l,f)l2  f  !  ip3 

1  “  I*  ,  ,,  '  irVir  4  • 

^dfjyi,1 


d-9) 


Since  s(t)  in  ( I -1 )  peaks  at  t  =  1 ,  we  expect  that 


A  cos(2irf1(t  -  1))  (1-10) 

ought  to  be  a  good  fit  to  s(t)  for  t  near  1.  In  fact,  plots  of  (1-1)  and 
(I -10)  for  A  =  exp(-.5)  show  very  good  agreement  for  .8  <  t  <  1.7. 


Thus,  smoothing  of  the  WOF  W$  in  (1-2)  by  means  of  WOF  Wy  in  (1-4), 
for  this  example,  results  in  a  very  meaningful  distribution  function. 


I -3/1 —4 
Reverse  Blank 
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APPENDIX  J.  DOUBLE  CONVOLUTION  OF  TWO  GAUSSIAN  FUNCTIONS 


By  means  of  the  double  integral  result 

rr  1212 

JJ  dx  dy  exp[-  2  By  +  yxy  +  yx  +  vy]  = 


[a0  -  y 


2  2  1 
,  Pu  +  ay  +  2yu«  I 

L  2(“B  -  Y*)  J 


(J-l) 


for  ay  >  0 ,  (3 r  >  0,  oyB^.  >  Yp »  it  is  readily  shown  that  the 
double  convolution  of  two  general  Gaussian  functions  is  given  by 


exp[-  2  a)<2  "  2  by2  ”  pXy^  ®  exPt“  2  c*2  ~  2  dy2 


2.  J  *  N2^  *  2H3x, 

■  jf72  "Pr  - 20 - 


-  Vcff  Xxy]  =  ( J-2) 


(J-3) 


for  a,b,c,d  >0,  \p\  <  1,  |x|  <  1,  where 


D  =  ab(l  -  P2)  +  cd(l  -  x2)  +  ad  +  be  -  2)fibccPp  x 


N1  =  ac[b(l  -  p2)  +  d{l  -  x2)] 


N?  =  bd[a(l  -  p2)  +  c(l  -  X2)] 


N»  =  Vabcd1  [VaF  x(l  -  P2)  +  P(  1  -  x2)] 


(J-4) 


J-l 
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Also,  a  useful  auxiliary  relation  is 


N1N2  -  N‘  =  0  abed  (1 


p2)(1 


x2) 


(3-5) 


Now  let 


p  =  sin(e),  where  -  -  <  0  <  j 


5$ 

fe 

I 


X  =  sin(0)#  where  -  |  c  0  <  | 


(3-6) 


Then  the  area  of  the  contour  ellipse  at  the  1/e  relative  level  of  the  first 
exp  in  ( J-2)  is 

2  ir 


A,  = 


1  Vab*  cos(e) 


(3-7) 


in  the  x,y  plane,  where  we  used  (D-l),  (0-19),  and  (D-20).  Similarly,  the 
area  of  the  second  exp  in  (J-2)  is 


A_  = 


2* 


2  Vcd*  cos(d) 


(3-8) 


The  sum  of  these  two  effective  areas  is 


A.|  +  A^  =  2w 


VaF  cos(e)  »  VoT  cos ( 0 ) 
Vabcd1  cos(e)  cos(0) 


(3-9) 


On  the  other  hand,  the  area  of  the  contour  ellipse  at  the  1/e  relative 
level  of  the  smoothed  exp  in  (J-3)  is 


A3  =  2ir 


3EL 


cos(e)  cos ( 0) 


(3-10) 


J-2 


IMUULMHIMHJWliUWU  IT.Wrj  EUUMVFV  SC  <L"  HFVKPMWnF.  V  f-TA'A' '  ?■.  -JWWW’JV  gg  UW  7TVW  WJiff 
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in  the  x,y  plane,  where  we  can  express  0  from  (J-4)  as 


0  =  ab  cos2(e)  +  cd  cos2(e)  +  ad  +  be  -  2)fabc<?  sin(e)sin(#)  = 

=  [VaP  cos(e)  +  VciT  cos(e)  ]2  +  J/aT  exp(  ie)  -  /UF  exp(  i#)J  2  .  ( J-l  l ) 

Comparison  of  the  square  root  of  ( J-l 1 )  with  the  numerator  of  (J-9)  reveals 
that 

A3  >  A]  +  A2  ,  (J-l 2) 

with  equality  occuring  if  and  only  if 


yaP  *  YET  and  e  =  e 


( J-13) 


That  is.  In  order  for  A3  =  A1  +  A2,  we  must  have 


~  *  ~  and  X  *  p  . 
c  a 


(J-14) 


Physically,  this  requirement  states  that  the  contour  ellipses  of  the  two  exp 
terms  in  ( J -2 )  must  have  the  same  ratio  of  major-to-minor  axes  and  they  must 
have  the  same  tilt.  If  either  condition  is  violated,  then  A3  >  A1  +  A^, 
the  exact  amount  depending  on  the  second  term  in  ( J-l 1 ) . 


EFFICIENT  CALCULATION  OF  GAUSSIAN  FUNCTION 


If  the  general  two-dimensional  Gaussian  function  in  ( J -2)  is  sampled  on 
an  equi-spaced  grid,  for  purposes  of  convolution,  it  will  be  necessary  to 
compute  the  quantity 


02(m,n)  =  expf-am2  -  bn2  -  cmn) 


(J-15) 


J-3 
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for  integers  -M  <  m  <  M,  -N  <  n  <  N.  The  following  efficient  procedure  is 
based  upon  the  general  method  given  in  [12]. 

We  observe  first  that 


02(-fli.-n)  -  Q2(m,n)  .  (J 

which  cuts  the  effort  by  one-half.  There  also  follows 


02(m.n)  =  Q2(m,n  -  1)  Q^m.n)  , 


(J 


where 

Ql (m, n)  =  exp[-b(2n  -  1)  -  cm]  =  Qi(m,n  -  1)  exp(-2b)  .  (J 

These  recurrences  can  be  started  with 

Q^m.O)  =  exp(b  -  cm)  , 

Q2(m,0)  «  exp( -are2)  =  Q2(-»,0)  .  (J 

Furthermore,  these  latter  two  quantities  are  available  through  the 
recurrence 


Q^m.O)  =  Q^m  -  1,0)  exp(-c) 
Q.(m  -  1,0)  »  QAm.O)  exp(+c) 


for  m  >  1  ,  ( J -20) 


J 


TR  8225 


Ql(O.O)  =  exp(b)  , 


(J-21) 


and  the  recurrence 


02(m,0)  =  Q2(m  -  1,0)  E(m) 


£ ( m)  =  E(m  -  1)  exp(-2a) 


for  m  >  1  , 


(J-22) 


02(0,0)  -  1,  E(0)  -  exp(a)  . 


( J -23 ) 


The  only  case  not  covered  by  the  above  recurrences  Is  for  m  =  0;  then 


02(0.n)  =  Q2(0,n  -  1)  F(n) 


F(n)  =  F(n  -  1)  exp( -2b) 


for  n  >  1  , 


(J-24) 


F(0)  =  exp(b) 


(J-25) 


A  program  for  the  evaluation  of  (J-15)  is  given  below.  Only  three 
exponentials,  in  lines  90-110,  need  to  be  evaluated.  Also,  the  only  storage 
required  is  for  the  final  quantity  Q2(m,n)  in  lines  60-70.  The  auxiliary 
variables  Q^m.n),  E(m),  F(n)  introduced  above  need  never  be  stored.  The 
check  on  accuracy  in  lines  390-470  would  be  discarded,  of  course,  in  any 
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practical  application;  It  is  appended  as  a  check  on  any  typographical  errors 
in  entering  the  program  Into  another  computer. 


10  ft*. 037  !  »xp<-8  »-'2  ~  B  n~2  -  C  m  n ) 

20  B*.051  !  for  -M<*(»'*M,  -N<*n<*N 

30  C  * • 044 

40  M*3 

50  N-7 

60  REDIM  Q2C-M:M, -MjH) 

70  DIM  Q2<30,30> 

80  DOUBLE  M,N,Ms,Ns  !  INTEGERS 

90  Ea-EXP<A> 

100  Eb*EXP<B> 

110  Ec-EXP<C> 

120  Q2<0,0>-1. 

130  E*Eb 

140  E2b*Eb*Eb 

150  FOR  Ns*l  TO  N 

160  E*E^E2b 

170  Q2<0,Ns>*Q2<0,Ns-l >#£ 

180  NEXT  Ns 

190  E-Ea 

200  E2a*Ea*E* 

210  Glpo*Ql»o*Eb 

220  FOR  Ms* 1  TO  M 

230  Qlp-Qlpo-Qlpo^Ec 

240  Q1m-Q1mo«Q1no*Ec 

230  E-E-E2a 

260  Q2<-Ms,0>«Q2<Ms,9>»Q2<Ms-l,0>*E 

270  FOR  Ns* 1  TO  N 

£30  Qlp*Qlp' E£b 

290  Q 1  m  “  Q 1  m '  E£b 

30O  Q2' Ms, Ns >*02 'Ms , Ns- 1  >*0ip 

310  Q2<-Ms,Ns>*Q2<-M*,Ns-l  >#Q1#» 

3£0  NEXT  Ns 

330  NEXT  Ms 

340  FOR  Ms«-M  TO  M 

330  FOR  Ns* 1  TO  N 

360  Q2<-Ms, -Ns>«Q2<Ms,Ns) 

370  NEXT  Ns 

330  NEXT  Ms 

390  Big»0.  !  MAXIMUM  ERROR  CHECK 

400  FOR  Ms*-M  TO  M 

410  FOR  Ns*-N  TO  N 

420  E*EXPOA*Ms*Ms-B*Ns*Ns-C*Ms*Ns> 

430  Error*E-Q2<Ms,Ns> 

440  B i g*MAX< B i g , ABSCError ) ) 

430  NEXT  Ns 

460  NEXT  Ms 

470  PRINT  Big 

480  END 
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